{
  "nbformat": 4,
  "nbformat_minor": 2,
  "metadata": {
    "colab": {
      "name": "scratchpad",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3.8.10 64-bit"
    },
    "language_info": {
      "name": "python",
      "version": "3.8.10",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "interpreter": {
      "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "source": [
        "import jax\n",
        "import jax.numpy as jnp\n",
        "import numpy as np\n",
        "jax.devices()"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[TpuDevice(id=0, process_index=0, coords=(0,0,0), core_on_chip=0),\n",
              " TpuDevice(id=1, process_index=0, coords=(0,0,0), core_on_chip=1),\n",
              " TpuDevice(id=2, process_index=0, coords=(1,0,0), core_on_chip=0),\n",
              " TpuDevice(id=3, process_index=0, coords=(1,0,0), core_on_chip=1),\n",
              " TpuDevice(id=4, process_index=0, coords=(0,1,0), core_on_chip=0),\n",
              " TpuDevice(id=5, process_index=0, coords=(0,1,0), core_on_chip=1),\n",
              " TpuDevice(id=6, process_index=0, coords=(1,1,0), core_on_chip=0),\n",
              " TpuDevice(id=7, process_index=0, coords=(1,1,0), core_on_chip=1)]"
            ]
          },
          "metadata": {},
          "execution_count": 1
        }
      ],
      "metadata": {
        "id": "lIYdn1woOS1n"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Sample MCMC chains in parallel\n",
        "\n",
        "https://github.com/blackjax-devs/blackjax/blob/main/notebooks/Introduction.ipynb\n"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Setup"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "source": [
        "pip install blackjax"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Defaulting to user installation because normal site-packages is not writeable\n",
            "Requirement already satisfied: blackjax in /home/kpmurphy/.local/lib/python3.8/site-packages (0.2.1)\n",
            "Note: you may need to restart the kernel to use updated packages.\n"
          ]
        }
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "source": [
        "import jax\n",
        "import jax.numpy as jnp\n",
        "import jax.scipy.stats as stats\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "import blackjax.hmc as hmc\n",
        "import blackjax.nuts as nuts\n",
        "import blackjax.stan_warmup as stan_warmup"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Model"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "source": [
        "loc, scale = 10, 20\n",
        "observed = np.random.normal(loc, scale, size=1_000)\n",
        "\n",
        "def potential_fn(loc, scale, observed=observed):\n",
        "    \"\"\"Univariate Normal\"\"\"\n",
        "    logpdf = stats.norm.logpdf(observed, loc, scale)\n",
        "    return -jnp.sum(logpdf)\n",
        "\n",
        "\n",
        "potential = lambda x: potential_fn(**x)\n",
        "\n",
        "initial_position = {\"loc\": 1.0, \"scale\": 2.0}\n",
        "initial_state = hmc.new_state(initial_position, potential)\n",
        "print(initial_state)\n",
        "\n",
        "inv_mass_matrix = np.array([0.5, 0.5])\n",
        "num_integration_steps = 60\n",
        "step_size = 1e-3\n",
        "\n",
        "\n",
        "\n"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "HMCState(position={'loc': 1.0, 'scale': 2.0}, potential_energy=DeviceArray(59465.977, dtype=float32), potential_energy_grad={'loc': DeviceArray(-2158.9343, dtype=float32, weak_type=True), 'scale': DeviceArray(-57353.906, dtype=float32, weak_type=True)})\n"
          ]
        }
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Single chain"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "source": [
        "hmc_kernel = hmc.kernel(potential, step_size, inv_mass_matrix, num_integration_steps)\n",
        "hmc_kernel = jax.jit(hmc_kernel)\n",
        "\n",
        "def inference_loop(rng_key, kernel, initial_state, num_samples):\n",
        "    def one_step(state, rng_key):\n",
        "        state, _ = kernel(rng_key, state)\n",
        "        return state, state\n",
        "\n",
        "    keys = jax.random.split(rng_key, num_samples)\n",
        "    _, states = jax.lax.scan(one_step, initial_state, keys)\n",
        "\n",
        "    return states"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "source": [
        "%%time\n",
        "rng_key = jax.random.PRNGKey(0)\n",
        "states = inference_loop(rng_key, hmc_kernel, initial_state, 10_000)\n",
        "\n",
        "loc_samples = states.position[\"loc\"].block_until_ready()\n",
        "scale_samples = states.position[\"scale\"]"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 992 ms, sys: 0 ns, total: 992 ms\n",
            "Wall time: 4.63 s\n"
          ]
        }
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "source": [
        "\n",
        "print(loc_samples.shape)\n",
        "\n",
        "fig, (ax, ax1) = plt.subplots(ncols=2, figsize=(15, 6))\n",
        "ax.plot(loc_samples)\n",
        "ax.set_xlabel(\"Samples\")\n",
        "ax.set_ylabel(\"loc\")\n",
        "\n",
        "ax1.plot(scale_samples)\n",
        "ax1.set_xlabel(\"Samples\")\n",
        "ax.set_ylabel(\"scale\")"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(10000,)\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0, 0.5, 'scale')"
            ]
          },
          "metadata": {},
          "execution_count": 7
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAAFzCAYAAACZyCAUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACCV0lEQVR4nO3dd5gT1foH8O/ZzgJL723pVepKEVApCoq99371qle93utV7Hr1p9h77+3au4AKgkgREKR3kKWXpffdTXJ+f2QmOZnMpCeT2Xw/z+NjNskmZ0OSmfec876vkFKCiIiIiIiInCnL7gEQERERERFR7BjUERERERERORiDOiIiIiIiIgdjUEdERERERORgDOqIiIiIiIgcjEEdERERERGRg+XYPYBI1K9fXxYXF9s9DCIiSrK5c+fukFI2sHscTsHjIxFR5gh1jHREUFdcXIw5c+bYPQwiIkoyIcQ6u8fgJDw+EhFljlDHSG6/JCIiIiIicjAGdURERERERA7GoI6IiIiIiMjBGNQRERERERE5GIM6IiIiIiIiB2NQR0RERERE5GAM6oiIiIiIiByMQR0REREREZGDMagjIiIiIiJyMAZ1REREREREDsagjoiIiIiIyMEY1BERxUBKidXb99s9DCIixztU4cLG3YfsHgaRozGoIyKKwYez1mP4079hTukuu4dCRORoXe77CYMem2z3MIgcjUEdEVEM5q3fDQAo3Rl6drnC5cGY8ctRtr88FcMiInKsI5Vuu4dA5FgM6oiIYiCl9/+7DoYO1m79bD5enbIGgx6blIJRERE5V3mlx+4hEDkWgzoiSmtjF27BroMVdg8jiNvjjeoeGbcc387fZHm/sQu3AADKXTxZISIK5TBX6ohixqCOiNLW9n1HcOP//sTfP5xr91CC7D7kDzRv+WS+fQMhIqoiGNQRxY5BHRGlLX11a9PuwzaPJNjUVTsiul/twtwkj4SIqGo4XMGgzql2HazAiq2sCG0nBnVElLaWbdkHANi0J/2Cuki1bVDD7iEQUZqQUmLdzoN2DyNt/eeLBXYPgWJ0yvNTMeLZ3+weRkZjUEdEaevaD9Jv26WuYc183+X+bepa3m/uOm+VzM5NipI+JnIGIUQLIcRkIcRSIcQSIcQt2vXnaj97hBAldo+TEu+7BZtx3BO/YlqEK/2ZZsnmfXYPgWK0ee8RAECl24Nr3puDf306394BZSAGdUSUNnYeKPcVIElnRyrdAeNsUacw7O+IZA6InMYF4N9Syi4A+gO4UQjRBcBiAGcB4HR3FbV4014AwJLNe20eSfrY7OCdGBTsnelrMXHZNnw1z7qAGCUHgzoiSgu7D1agz8MT8cRPKwAgZEVJu53x0nTsVCpyfj53I/5v7FK8O32t5e8s3cIZaPKSUm6RUv6pXd4PYBmAZlLKZVLKFfaOjpIpN9t72lXpZjVcXTpWN6bY7TpYafcQMhaDOkqadTsP4vSXpmPPIX5hU3h7D3sPBK9OWQMpJRZs2Gt6u91cbg+WmySDvzF1LR74fmnAdVKm/6oj2UsIUQygF4BZUfzOtUKIOUKIOWVlZUkbGyWe0Jbs+dXgt6bsgN1DoATavu+I3UPIWAzqKGlemLQaCzbswc9LtoW838y/dqLk4QnweCRWbtvvy0GizJIl/BsUy/aXw+UJnMleniYrXUei6Dc3ffXOgJ958kIqIUQNAF8C+KeUMuI3uJTydSlliZSypEGDBskbICVchfb9UXagHG9O/cu3HTOTsSVM1cJtl/bJsXsAVHXpM5HCIpno2McnY/2uQ76f7/5mET6evQEAUDpmVLKHl7F2H6zAjgPlaN+opt1DCeBRpq4PVbjRrWmtgNufnrASn143INXD8jn31RkorlcdYxdtifh3LnkrcPHlpUmr8fT5PRM8MnIiIUQuvAHdR1LKr+weD6XGG1O9W7Tf/32d7zoe7/yO7+ifpCh3uSEgkJfD9QeiSPCTQkmjbz3Lsojq1IAOgC+go+Q67aVpOOGZ9KvDUKHkmJS7PEEH8llrd2H7/iP4fsHmVA8NAPBH6W58PncjDsXRR4kzmAQAQggB4C0Ay6SUT9s9HqJ04fZI3PnVQpz36u/ofO+POGbMJLuHRBHKMjnVe3XKmtQPJIMxqKOk0RtH//vz6PvOMIk8eTbsSs9KYz8s9K+AzVq7E+Uub/B0YpdGvutPf3E6bvp4Hsr2l6d8fGb6FtfFW5dHX3n+t5VleOO3v5IwInKIgQAuBTBUCDFf++9kIcSZQoiNAAYAGCuE+MneYVIqMPfWb9Oew/h49gbMLt0FjwR2HEiP73oKr7h+dQBAjhLdjRm/3K7hZCQGdZQ0akWrA+UuPPj9EhyOcJXjrzI2Z022dGsdoJ7YtKxb6JsUuHtUZ9/1W7Q+OIcqXCkd2zKLfL6nzuuB+jXyA64zKwzUpkH1gJ8ve3s2/m/cssQNkBxFSjlNSimklN2llD21/8ZJKb+WUjaXUuZLKRtJKUfYPVZKjqa1CnyXfwqTd15VPfj9EhSPHhtwHY/9zlU9z5vR5Uqzc4tMwqCOkkZdbXttyhq8M70U7/9eGtHvnvRc+m0PrGoOV8a+jTDRyvaXo3SnfzvuwXI3yiu975+aBblB949nC2QsTnpuatB1X15/DFrULUR+buDX6AqTypjvX9UXAHDT0HbJGSAROcIp3ZsA8DdqBoAv5m60azi2emd6qd1DoCj8uHgrvgzxXuUOK/uxUAolTYOa/hUMveKXO8JtJh4JFI8eiyUPjkD1fL5Nk+FwhRs10uC13br3CPo/+kvAdQcrXNimlUU2G2OqgzozfVrVAQDk52QHXP/ETyuwcNNe3D6io++6omrewLRWteAAlYgyR7Pa1YKum7gsM1fqyFn+/uFcAMDZfZqb3l4RRWVoSg6u1FHSNCrybi/p0aI2XtPyh6JNHdiyNz3zv6qCycu32z0EAAgK6ABg3vo9eHOat0pcbnZw9vXhCjf+KN2Fk56bin1HUt+/7qlze/gu5xsKusxZtxsVLg8eHuvfXpmtFQt65dfApPFnJqzEm1OZW0dU1VjlyaXbtvd08h9lIoyc5UC5C3/tMN86y/d86jCoq0ImLt2G4tFj0yYQ0qsZlivb/PQT8Eg/5JVufhkk0urt/j5pt3+5EFNWpmfj4o9nr/ddFibVU1dt349zX/0dy7bsQ/cHfk55oYGTjmrsu1yQmx3inl7ZWuL4zoMVAXmlz/2yCg+PXZY2jdWJKH5b9h5G6zvHmW5VY76RXxMlr7BmQQ6Wm2xdJ2d4afJqy9u4Ep06DOoc6oHvluDJn1YEXPfJH94T4fNfm2nHkILoOVHqF/X6nYfw9M8rsN+wunLHyE6mjxFp8DendBdL50Zg+NNTAn6+/O3ZNo0kPsZ2GMkMijbtCZ4kycv2f3UaV+qMTu/ZNKCtx8KNe4Luc+ZL02MfIBGllVXbvJNnX80LDuo27j4UdB3g7duaafq2ruu7/K8TOli2qzmSRvnfZJ47Z1YgTPfUzyssb6PEYlDnUO/OKMWLFjMjxhNeu1SYfPDHL96K5yetxouTAsdeq1ouLu7XMubnOufV31k6N0EqXB5c8Prv+HP97qQ/15qyAwE/P3b2URH9njHYX7hxb8LGZLT7YPDBKieKoG7uut0BJZ7Pfz140sVq2woROY8enExfvTPotonLzLe9p8txO5XU4MAVYlcOVzfTi9mqXKg+w9x+mToM6hxO3T6Xbq1uKlzWs2t6vpQuJ0vgnlFdgu5XHmXiLfv9WNsZYb+fdTsPYuZfu3DWyzOSOp6v523EsKcCVw7PPzowsO/fpi7MGFtjvPZb8lZps0y2f6rUAM/M1YNaIytLoHvzWijI5VcuUVW3aJP5JNMTP3HiUVXh8h+v9cJYZr76MzOrg6arZyeuCnufVvUKfZeLWCAsZXiG4UDqMveUFf6gzqpE/ZzSXViwYU+yhxWg3OWOqvdOYX42quVlo3TMqMDHiXLbRbRBYKZwuT34dr751hYjYzXHZLn10/BN6bdrTcaHd24UcP3UVTsCfjabEU+Ug1pPvFgrVx7Ttj4A72rikUq+P4mqulxlokedaHxpMlMEdFLKgFwrl0firF7NTO9737dLUjUssvDpH+vD30lz5cBiTPzXcb6fmTOeOgzqHOinJVt9l10e/0nijDX+E9tv52/yXT7n1d9xeopzdn5cvDX8nRTHd2xoev2TUe7F3nWwAgfLU9uY2gna3T0e//1haUT3Vd9T5SFWW+OVnRW4AvZ/Z3YLuo/eiLZdwxoB128NMaubaPd+sxgAcOOQtvjHkHb46oZjwv7O+1f1RZv63objYRbyiKiKOe/oFr7Lre8cF9XvSimxy2TLd1VyqMKFuesCt/fn52bhPyNZ/TJd3fHlooju9+HV/XD/qV2Rm52Ftg28x8C/yg5W+fd0uuDphgOpOwyPLq4LKWVQWfdIV2WSJdo91Fb90v5cvyeqxzlmzCR0vf+ngKA2080Lkxs3cem2gEBYzV8wHngTyfgeOambtynvO1ceHXTfk5Vqk6q7T+6c+IEZ6IV+3B7gthEd0btlHcv7zrprGO4Y2QmD29fHCxf1whk9m6K4nvfAVtIq8Pcm/fs4s4cgIofzGL7bIk0LcLk9+HzORvR+aAKWbt6XjKGlhVNfmIZzXv3d9/NtJ3bALcPaBxSgImca1L6+7/Ibl5X4LpsVCKPE4yfIgfYoS9nV87PR+s5x6P7Az2hc5C8PrOfuPDpuWdDvA0DpjoO47oM5SenxtfdwJf71WfitdQBQOmZU0JbLe0bFf6I+btGWuB+jqvhsTuh8hGven4Ou9/8EADhY7sKJz/zmu00gdD5ZIhXmebd96kGQqnvz2njojOCVvNb1/fc1nkgl2mDlYGV01cDWeOOyEjQqKsD1x7eFEAJdm9bCsxf08uXczTEEyOrYAeaDElUVhww5v8u2RFaqv93d4/H5XG/BiZXbqm55/zVlgYWh/jG0PQrzcpAbpugUpSdjNXNd09rVfJc/nLku5sffe6gSD/2wlFVQI8BPkAN1bFzTd1lNNK5dmKtc791Cpzf9BgJPGseMX46flmzD1JWBuUnxcnskejz4c1yP0aZBdCe7Zvu1o8nnyzQ3DW1ner2U0rYKbBf3a+mrImmVu3Zp/1ZBEwDV8vz5f2bVVhPJuAVUdd+pXXBCl0aWt5sRQuC2Ezv4fv59TfLyAokodYwl36/7cE7Ev/tHqXfyx5OBkzxcqXMmtf+tKjC3NPbHv/yd2Xhr2lp8PY87sMLhJ8iBPlVKx6r5T+qe5YnLtgf11/pcWbH5UcvL23M4cfucN+4+hLZ3BeYPdFIC0Egd16FhwGrdFyYNXHUd7hkfdxBZ9QV+m15xTLHpvc5+ZUZQM/JJy5MXHI86yrvd8qlze+D/zjzK12S8MC/yQi0D2/lXz5JRJEf9TOUm8ITjzpO8fRnVxuofxDGTSUTpwzhBs2FXcK/LcDKxjH8iv2MpdSYsNT9PUPPmz+rdPObHn5/iQn9Oxk+QA3VpWuS7rM4I6pUCdQPHTAr4+fYvFwY91ncJzL373GSb372ndMHbV3j3Vc+9ZzgWPzgi7ONkZwlcM7gNTtRWPooN29QmL9+O4tFj8ePirb4VSbKmruYCQJ3CPPxw06Cg+/25fk9Qr7+dSUhu3nOoAk/8tBxjtS2yZ/cJ/LJX+74ZtyiGEm2l1Ehc8Lo/78NY2CVa1wxq7busb0s5V/nbxy/eilL2qyNyvN//Cl51V1vKTPzXsWEfI9nbydOR+h17fkkL0+MU2UsIoJmyrRIA+repF/b38hKwtTYnzmNwJmBQ50BNavlz5776M3A5umXdQuPdA7g9EsOe+tX3844Ie5dFoluzWkHXHdO2HoZ2aoTSMaNQr0Y+auTnYMXDIyN6vKu1k2Bj4Hblu38AAP7+4dw4R5wZvjT0+MnKEujWrBZ+ve34sL/bu2UdrNi6H/+bFXk543B6/ndCyNLe6urVF38fEPbxnjinO4DkrNSt3Ga+rSQWA5WcvFVavkzDogJMVv4djn/y14Q9X7z2HqrE9hRWGSWqytSWQ5H07crElTrVY+d0R8t6oc9nKDVcyuJBfk5WwA6xcPR0D1cC0iMSERhWdXyFHEj9sp+6akfAG32zYculcXVh2uodAUnKxoTleJg1VhYmjZsj7YOm/11VOWE82UJt6TOugJrZe7gSI579DXd9vQjjE1R8pnoE2yufPLcHvrlxIOrVyA+6beXDJwX8XJDrfbx0T6LOzfJ/PoYqfffyDQeqvYfSo6dPj//+jL6P/IIpK8vwdJStRYhCWbhxD6avTmw+d7q74PWZALytUdTvAitVLafuYLkLhyqC2w0Zc7zfu6ovfvznYABAUUEuBrWrj3rV81IyRjJ3r9InsMLlwbZ9gYsBoXZMnd6zKYDETFLkRPC5yXR8hRzIOOOhxm3GD46xbPzlb8+2fNxO947H87+sinlcxg/21NuHxPxYALBHO7l91LAl0Mqf956Afq3rxvWcTudye7B4017fwVPvsRarJ37yn8xf/9GfcT2WTn1HXtq/lel9zunTHD1b1Da9zThb5w/qErtSp37O+rSybmMQqZxs/wdVLQhgPNStLgs/ibFw4x489fOKqFuHREp93Mvfno3nJ61O6Ko+ZbbTXpyOi9+cZfcwUmrjbu+E67It+yNacXC5q1ZQ1/X+n9DnoYlB1181sHXAz8d1aIBOjf0pJsX1uVpnt8/n+Os46IcG9Rjxy/LtAICaJq2psrVALJrVPdVnf/if2zgBSsGS9goJId4WQmwXQixWrqsrhJgghFil/T/+M6UMpM6aAIEns9VyIy8yoXtvRin+KN2FI5UePD1hZczjUoO6G45vixYhtoI2r1MNF/ZtGfLx9CIYfz+2TUTPX7d6Hj69zr9dz1iBLBMMeepXnPLCNHS57ye8O31twG2/3nY8lj8UuPX1pYt6p3J4AALLfRdVM+9PGA19hTiRDcnLXW68rbx+p3ZvEvdj5ipBnbqCbpyFPlgefsXxtBen44VJq/HVn6HbVcTqHcN7BwBKHg4+ISOi6ORkCVTPz8H9p3bBe1f1RYOawbsRAGDRpr0pHlny6PmBh012U2Rnh86TysnKSnplYwrNbJWt3OX/t/x4tjc945Pr+gfdT8+Dq4xxkkKtBVHVVq+TIZlh77sAjMlTowH8IqVsD+AX7WeKQrjy/ocr3fjomn5RPeb93y3BuUoj0FipX7zhmlZPu2MoHj3rqJD3ycvJQk6WQHmMX+jt7x5vut2jqvpuweaAKmsPfL804Pbi+tV9q1q6UQkIVuLRuFa18HcyceOQtr7CI/oK0t/en4PDFYnZgnnpW7PxyDj/CvHlFhVDo6FuRa6jtB8x/ptEMxkRaf+raD081ry/JRHF57rjvJOUVw5sjeM6NMDT5/VAK5PcsapUvn3ZVn8jdeMxOVzxi7ycrKgnaP8o3YVtzAdOCn3hoNItccsn81A8eqzvNrOWFHpF00TsKmFQF17Sgjop5W8AdhmuPh3Ae9rl9wCckaznr6oiKQahlnmPxVvT1sZUeUsd27qdiel35vJIvDbF32vPrCedkbqd48AR5wd1Uko88N0SLNuyL+T9JpqUFT6vJHwZ4dcu7RPz2OIVyfjM/GdEJ9xzShcAwDFt/e/3q9/7IyHjmr028KvLLDc0WupKdsOigoDbPr22P07r4c09MDYuNlIndsxOBs3uHyrn4eVfV+PtGD/zRBSsRn4OaphsRQOAo5rVDvh5cPsGmPKf+FIV0p26867LfT8F3BauqnButoh6K+q5r/6Ok56bGtXvUGTO0So2V7o9+NZQPd3sX0n/97XqZRcNLtiGl+oNqo2klHq1ha0ALLv1CiGuFULMEULMKSsrs7pbxlGLQbRtEHm591O1E8ZIPPTDUkxYFn1/sseU3LcTu0bXiDlSq7cHr0wM79wIM+8c5vu5eR1l9acKVMA959Xf8e6MUpwXZjXV2JcQAD4zaTNhdKJJ0+z/nt418gFGSc0nibRoTiiNlABpRho38NZnm3s0D64S269NPfxnREcACLvaqPbAWrBxT9jnfeKnFehwz3jTQjIVLg8e/3EF/vvDUt92U2OvQtUGm5rTEzmJy+NBNYuCUNFU8Fu6OfREnlOEWmEJV/wiJysLLo+MaNLp3m8W44HvvOkpu5LQjkf16wpvayX1nGTfkUq88MuqpOU6p4PGWvV1s0C7qCC4qqse1L01LXhLfyhSStzzzaKA69xcqQvLtqxD6Z1utvwXklK+LqUskVKWNGjQIIUjS29q/lxNkw+QlUfO7BbV8+yLYEUM8O6V109C1Z5midiuBngLVLRrWMP389mvBAc2d57cyfdFA9hf9jbcFtlo6VtZ95eHXnUMteV19t3DLG8zW4W6bEAx+hYHF53ZmYBiGeHabqSbly9OTN5hcT3vJMxF/cxzSfWTQLO8E9XkFdt9l40tTcx8pLWjMAsWlyvbon5esg3d7v8pZDGl+76Nr/AOkSrdK9bG4uaP5+FIpcdX5XdIx8jOX7664Zig605+vmqsNhnb6qjCrdTpx/PKCAptfDBzHd6dURrV2GI1ftFWAMDwp3/zXffQ90vx1ISVmLR8u9Wv2cbl9mB3AgJdPQfUbEuseh6mi7W1nNsj8eHMwFZK3E0SXqrPfrcJIZoAgPb/9Hvnpzn9IDjmrKPw5uUlQbc/eW4PAMDEfx0XcH00ASAA/OeL4EblZsb8uByd7/sx4OD8/T8GoW2DGiF+K3INa+Zj0+7DpqtQgLehs/G51ApJMRZcitnpL05D6zvHpd1+/oY1g79swzELjj+fG39hjtXbD6BLk6Kg92i6WLczsM3HCSYrmbFoUbcQyx8aifOPNg/q9Nf7/u+W4GCIAP5glHmD+kH1iCv49057cbrv8uzSXTgQZuIgXMBJFI13ppfaPYSE+26Bd0uaXhhCzZk9JUQOc9MY84udwOp4GC6gA/wFpmKpBnr7Fwui/p1IfWESqO7WKnanS+7XlJVlWL19PzbuPoRHxy9Hr4cmxPV4X14/wJc3F2meY+1CfyGwcMcXlVlhlUzv3RiJVAd13wG4XLt8OYBvU/z8ae1guQtjxi8Puf1KPzErqpaL+oYeXpf0b4kzezUDAFTPD7+tLRG9Xz7TSt0eLHfh8gHe8vRHmWwvi1V+ThYOV7oxcMwk09vfNFnSV4ORVAdXCzZ6K5bNWx+6UEy6uveULrhPy1UzC+rinVnX39tLt+wLWIGN14sX9fJdnvlXfFsw15QF7v3PNUn+jpWxKIpKTTKfv2GP5f0ORXFgBPwnGgMeNf8MRfVYB72PdaTSjb/KEteYnTLT/A3+78kZa3bgkjdnVZmqxeoWvPNLWgAABrStZ3n/qtyCK8fiO7Qggl01+vbMWN4XkaQfJNJELW0lN0xFz1TweCQuf3s2hj/9GwY9Ntm3/TGWraGdGtfEiV0aoU+rur7WPMYAyyqHVFUZQU0InVkeOFfqwktmS4OPAfwOoKMQYqMQ4moAYwCcIIRYBWC49jNpXpuyBq9OWeMLlMzo2y+Njb5P7NIID59xlG/mqzAv9Ads4r+OwyuXWBfIiKRBNODvJXf6S9MxNwmBTGkMBVfUE+dUb4PQqxpalamOhscjccZL08PfMQEu7d8K/dvUxdWDWuMqraqkWW7VsxNj72MIIGmlqU/p7s8ZjTeQt2uSVQ3qQgV/3xiS03cdrAhbbTYWl/T3ryjq7Ufq1/ROBA0cMwlDn5oSckWRyIx6YrZKKZ7w4HdLMW31Duw8kNxcqFQ5oUsjnNWrGe47tQv+fWIHnNStsa8YkhmzfCTAm8e6ctt+3Pb5Asfmao1duMX0+i+uD95yapSrBX7p1tYg1L9lIicCY/HJ7PWYvmaH6W2himZZqXR7fP8O+t9Wtj8wFSOSOgrRvHt3HPQ/fv0a3uMOc+rCS2b1ywullE2klLlSyuZSyreklDullMOklO2llMOllMbqmBlh0ca92LjbG6ws3bzPV2lvh7bfOdSMlL7SUWAoMDGia+OAn2tV8x8g9AIMf9w9HABwRs+maNewRsi9zkM6NYzkT/HZuPswFm9KfFK3WvFxxmr/l9T9p3ax/B31hLhzk5oJH5Nqx4FyvDejFNv3HcFHs9ahRoE3mP5+gflBLBqHKt0BKzaRVDqM1UNndMMn1w4IuE49gbj75M6+y/HMlunvbX01MJHO0lap1fd+tLbtO4L1SjEQtfVAsmUpH8hQZb7rGlbYz311Bs5+ZUZEq6jRVCDbsucIxpx1FDo3KfLl5E5f7V0F1fNnN+xm4RSKjpob9VeZf6vznsPe91S6bF2LV+3CXDx9fk80qVUNDYsK8MolfUKmQVhN5Ax+fDJOfOY3fDF3Y9AugnQjpUTx6LEBZe6tfHn9MejcpCjs/fLi2H6ZDJ/+sR7Fo8cGBTVqLn244i/JNvqrRbj0LfPc6Gmrd6B49Fi8/Otqy99/dNwyFI8e6zvWryk76Ks5px+bLn5zVsDvtKgT/vzEFSYwL91xEFv2elNt9AUDwF9QzamTGqlUhRf8U2PdzoNRV4Q79cVpGPTYZADeROjzXgss/vHw2GWWb159+2WBYSXt7D7WpeFvHNIOgHf1qHTMKDx7gXerWpbJiWOh9riRnPztPRRZMZV4qG0SLlK+RM7o2Qy52QLf3jgw6HdWbfNXo3px8uqEFy5R3frpfNz/3RL0feQX3P31Yl+fuHdnlGLL3sMht9GFY9yGG8mXJgC0qV8djYuiz6ELpU9xHd/lg3H0/tODOqvKcPG4WFtZirX9QIXLg36P/IIHlf5+3/1jUELGFq1QrTv0qm569ds12onx6C/D58EOf3oKfltZFtFn4pfl23FB35YYf8tgCCGQl52FmoYtNk/9vDLs4xCprI5tAul18h6vSPLFjH64aRBGHdUE428ZbHr7nhQcc+MRScslXZ9WdcLfCZFvv0zmcV41+itvRcZpqwNXwtT3tZ3bL8NNus7S0hMe/3GF5X1e+83bRurlX1f7CnP9oK22Wh1ebxraLuzYxi0KPdl9/JO/+tIECpVzBP18oapM+CQTg7o4HffErxj8+OS4H6d49Fj8b5a/0o9VQml5pflKXSwKTU6ss7RP7PKt4Zsa9/jvz3GPIZymJtWUAKBO9Tys+r+T0aNF7aDb1G0aizft8530JsPUVeZbHI7v2AADHp0U1/ZJY5PWtTsOov3d43DLJ/NCHuC+uuEYtKgbf9J916b+WdQ1SpAfy/YNXaXL+6WcjO0p+mxeeYx5f0/9HHiQ++/pXdHCpkqdxrGYMb6v1W2ZG3YdgscjTUuiX/b2bHwRQcEb40lphduD/eUunK68p6M5iSMCwhc7SLdtdpF6c+pf+Nv7c3w/h2uqbaZbs1p46eLelitYd0QwcWMnq/60er59LPRtf2aFM1RWO4W+nZ/YJu7HdQiuZrp5z2GMX7w1oc8Tq3CfL7UOgVmLKNX3C7bgyncCe79aTTha5U2qokmnUc9x9Ibn9327JOLfTxf7j1Tirq8XpSxVgUFdgsQySxTqd4wn9DqrnLpYdGocfOB47oKecT9uIsXSnuDY9oFfuomc3bnq3T/Q/5Ffwt7v1xXx91Y87olfA37etOcwKt0S387fjM9NEsCr52XjmkGtUbswD38b3MZ3/Q83xbbatEQJCPQCPADQ5+GJuCPC6qhG+glbMmYy9fdKrEH8R8qkCgD0bhnZTHIyGLdYmtHz3IzWlB3A4Mcn45Upa3Dqi9NM77Npz+Gg7UNG+RafvQXK6nO47TRERm6Lk3N9BSCeSSM7PTx2GSYs9fd3jeQkN5ROjYNTBzZbVIFOBw//sBQjnvWX97/+w7m+E/P3fl8X8+Pq2y9DTWS6PdLyu+6WT+YnfdveMWMm4aaP5/l+fuB7+4KPRZv2RnzfL+aGDnhXbAsO+nbsjz7ndXhnb75dNO0m1H9vJ7c+eernlfjfrPV4O8o+fbFiUJcg+oz14k178dLk1RH183p1yl+Wtx0sD34TezwS//x0PoDQhRR0424ejPeu6hv2fqphneMr357oHnGxbGFpXT+wKXu8Qd3va3b6vmAmLd+OrVEW4khGxaZyk/L0Lo9EtnYA1HP7AO/sbyz05GQg+ATl0xDFfADrLVb665iXlJU672M+9uPyqCdZtu49ErQ6vnG3fSdQQyPIae3dsrbp9Xru6e9rdlr+O2QJgaP/b2LIxzfrOWQUSU4MkcrYb0z/rOrf9FWl+mWk2+WtmBXbSuecImMl6vGLt+K8136P+/gXyfbLA0dCr4KYHS9jtTuCLbCLN+3Digh2OyXD2a/MiPi+7/9eGvXjq+cWkRoRQREVtZqylBIVLvOCSk6jB7JLNu+LqqVDrBjUJcjuQ97Zi1NemIYnflqBPg9PxCu/rgm6n7oEO8OiOhEA/LayDGt3BK44PDvRn78SSVDXpWmR6VaBRLA6af4gyiAyHLOE42fP7xnyd6ob8n7Uhu3RmrtuNy58Y2bAax+tWLYTvTTZOokZ8G+T1U1avg3lLo9vy8+ANvXwyJlHYcmDI6J+bt2P/zzW8rZQK0lfzN2ItneNM+0t6AvqktAgXn3MaL88Hxq7NOi6Do0S13IhWlbbGtW/y/g+B7z5dvdqW1RCbaFWV0P0Qko6fcvzqyGq4+oSnbtJVZ8xMHlhkve7Ts+FrSpB3clHNQ5/pyg5rU/XvPV7sO9IfHmAkWy/DHeMrXQl7nVbEGGe/PKtiS8cl2iHTNpnmZ23AkCjIu8kQyJbEanU87RKt6wy3wO6H5dsxYWvz0z68zCoS5CJS7dhv+HL67EflwfdT10at8rHAoD//rAUQ578NeC6sUqSqb798s3LSvB/WmW6eP3vmn4AgH6t64a9r9XMiVmOWzzOO7pF0HWRrCC8pTRmj2fpXu83ZwywoxFL3tETP4XPqdJt3H0IV73rzeXQvxiFELioX0vTE/9IGfsgqqxWUC99axZu+9zb8HVt2UFs3XskYMtQpW/7ZeK/etRAN1SvRzNmeXh1CuPv4xit207sAMD6JGXPIf/WF7NiM+rs8A7DbgG1yMmLyqSBmlt7Vq9m+OCafvj3CR3Q3nDwHtAmuMfW/41bZjpOIivGQij6NmD9PT9xWeRtaHYeKMdb09amrEhGpHq0qB1zwSansjpvCFX0KRK5WeGD/XABwOdzQ+8sSQYn/Puru3F0ZuetAPDuld4J+6OLw58fGp0aov2D/tk9+fmpvut2HCh37DbsUKLZGhsrBnVxULcVNK5VDWPGm38YVLEkT+tO7+nPa9ILpQzv0ggX94s9CVmlz4jpe/lDJXYaZ3jeuKwE0+4YEtEKYjSuGliMFQ+PRF/li6RhBD3g1G2k8QR1D4/1n7SqWzhcbk/EXzrRFO74cfHWoJnAOfcMx0hDywo130mvpArA12A0Ub68/hhMvu14AMCQjv5VX71twL4jlXh1yhrfF7M6UZElgP6P/oJjlMbxFUkslKKuAOiB9I4D5SgePTbs7KrZjGWODRXMbjjeW0FsikVOph64/uuEDkGrtQCwcOMe09+7/vi2vs+3UW52FsbfMhjjbh6Mp8/vibYNauCmYe2DTkq6NDXfavnAd0vwwi/x9S+kzGFcqftgpjffSg/uXp1ivlJg5t+fL8BDPyxF6zvH4V+fzU/YGOPVrHb8K9hOCApUbRpUN71+v8nWyHC7bVT+lbrg4+2hChfcHmlaMVWdYFaP4/FYH0WhDyf86zWtHXlBtaa1Yi++ZnZeeKDchRs+movWd47Doo2Bwc4pL0wLmti8uJ83h9zJ+XWpwKAuDmpzxCOV7qBCC0BwPtXhKN+QaiCh9swya0cQi8+uG4BL+3uDQr0P2h+l3tWpULMKxjy1E7o0QvM4cwjMCCGQn5ONx87p7ruuTYPIlv8/0lYeE/ElMG7RVsxfv8f382u//YUO94yP6Hej2f759w/nBlQXBLwrZj8uCaysVe7yJDRPwEqfVnV8OYq3DO/gu75/m7qocHnQ/YGfMWb8csxaG9xyUs//VPm3Xyb+kKfu9df/zU99wZs8f6vJWFQz1uwMus6OBrL659rs9QT82yab16mG9TuDV48ftZhYOqpZrZB5jJ2bFFkGbbqL+pkXZnl3RimemsDWBhQZYwufeKj5TV/9mdgqh+Gs23nQsqCZSMAp/TFtg1fGAQScAKfTCqXH4jCnD3FU9ya+685QCm+Fo38PGwM3KSW63PcT2t41znRnQ6EhkDDupIrWjgPlOPaJyCudO6GH56ndrVfQjNRMmKOVFkex6nb/Txi3yHteYyxys+tgBbbsDaxdoJ+HVLVtmYnGoC4Oz0zwz06vMqkSBARvo4o2v+uQUjClIgn9e/q2rouHzuiG0jGj0LCmd3bxH1q/kaIQjVKr5/lPoM2W8BOtuF4hrj++rWX/HjP6LFS4QHrv4Urc9vkC02qAHRv5K5Cdr+yHjmZ75O9/WW+zjdX93y3BTf+bF/6OCdSzRW2sffRkAMCHM9ej6/0/+m6bUxochGxXXk+9P40elCQjYFLfr/rnTD8w7D4UfcWueFbV45GXk4UzepofbN+Z7l2J3bL3CEqi2AZzUrfGyLUIpCM9+bDr9aiKhBAthBCThRBLhRBLhBC3aNfXFUJMEEKs0v5vXwnWJNh/pNL3vXDNoNY2jyY+xz3xK655b475jQn4qFx3bBvUM8ldPvXFadit9ars/sDPeDRNtkBb9R6tcHuPv4Pb1Y/pcfXvHeO5lPqz2Yn+akOj9gUbYt/65nJ7fN+9qpO6WedNPv7jCqwzmXhLJ8aiRaGoKRdDtEJeseRUW02EGC1WFhXysrN85wxVcVtmIjGoi4P6RfL8JPPCFsYvIrMAo3OTInx1wzEAgDFnHRVwm/rFpD/fJf3NZ8wTRc/XC5V8rP7tD5zWNanjAbwrdneM7BRVtb3q2kHGrJKo6sHvl+CLuRtNqwFabSkxY/UFf8eXiyL6/feUcr9tI3jen5Xy2brbR3Y0uWfiqFuC1MT1gtzskFXObvjoT5z58nRco/Vx2rYvfHXYWPz7BO9qonEVM9L6Auo211gqryZChcsT0HNOpZcG/3Dmuqg+C0IIyxYG+REG2GbbPVU/pkmfJodwAfi3lLILgP4AbhRCdAEwGsAvUsr2AH7Rfq4y1HYjxq1fR8VYpdcO+nedvsJvLHp2Ypf4qkgD3s9sW2VXSk1lJ0KvhyZASon95S5fo2i7WVXm3H3Qu0KWo23znvKf46N6XL0A1tRVgVvSS3f4J6PU85EzejbFz7ceG/R9d8lbs6J6XtUdXy7CS5ODtwVbrabqznk1cavS0fr2xoGYfdcwjL9lsGXf2mgKyKjf//oqaFG16HP2jbUirKjN3Uef1MkX1IXrV5hOtu6NrlJ6IjCoi0MvpaR439Z1fbMW6gn5dEMxFLOtgJ0a10TvlnVQOmYUzjcUBnl3eimKR49F6Y6DqNRmKB44NblBVCRNnE95wb9cblahMh0UaEFduO2XoXIHo5kVevLcHr7L0VYdnb12F+7/zt/bxthr7dpj2xh/xUevSgX4c7Ls0OaucSFvn6dsX03W1tF+WjEP44p4pNuUhndphLN6e7cGpVNOy5a9hwO2D0XTpuPCvt7vFKtdAu2V1ehQwlXei6dCbKaRUm6RUv6pXd4PYBmAZgBOB/Cedrf3AJxhywCTRF3tNU6a1NAK+QxuH/mKziabtripKxxSSlz0RmDAoOa/x0M9tBqPKepJbzr4ZVnwJCMAPPuL93shJ0ugc5MitKoX+UQp4N/V8eHMwPSW7fv9J8zXf/in7/KtJ3RAhwi/0yIxefl2fPlncF9YwH+uZCVcL9Bk6tGiNhoWFaBzkyJMvX2ob5eNKpp0IHX7fqH2WV257YBlP1MrkU7oqq/dJf1b+YJ7p2y/3Hu4Ev0fDd/TONHS82zcIe7+erHv8oA29TC4fX00KsoPaNRcz1BBUK+mqFKDQOOJpF7xcsLSbahwe5Al4m9qGo7+4bHqx2I8QU5GI+lE0L9sQlWfrHB5sFqp5GlspvzL8sgqsQ3v3CigiuAbl5XgkTOPCvEbgazyTPQA466TO5vevmXv4aStekUj2kT0UUc1CX+nGOirzMZAPtIc1NN7NsUT5/TA8odGJnxs8Rjw6CSMen6ar4rY1zcMBAD89cjJWPFw8FjV/JX/nm5dHffHfw4OWZlM1aJO6ET5UC0UyJoQohhALwCzADSSUupljrcCiH/JJ42ohzd1S7SUEiu1FIZoerHtOBD9tupEUPO7rPJfE2HmX97Hrl8jD8ZXRS2mlA65dZstViUWb/KW9o+2v6tOnQhQ/061uq/ePue1S/tEHTSGc+W7f1jepgbWxtYw6UY9t7xyYDGAyIsSXTmwOOAYqk6E/3nvCVgcQeukywfEXtAvLyfLd54ZSzXxZKpweUwXBqwWE0K1MksEBnUJ8twvq1Dh9qAgNxsNlX3GSzYH7uP+eHZwad1rBluvwujcUqLC7UlJ8QZ9yfiOLxeaj8Vw0J1tkk+VDvSZpXKXByu27se93ywO2CLo8Uh0uGd8wKqYuuVUnQkM9zyvX9rH96U5oE095OVk+U6s+0bQIsLK0+f1DHm7GtBF8sWaCnrhnVCStQqmV9kqd3nw8Wz/zG5JK/PUJLdHYu46//s3NzsL2Vki4VVco3HLsPYAgk/U1u865Mtd1LeuZWV5Cwl9/49BAfetoeS8hvrO6NQ48i2cOdlZKB0zCmNvHoTlD43EP4e3j/h3yZwQogaALwH8U0oZ0NhKet8ApmfrQohrhRBzhBBzysrMK6WmIzUPvdLtQXNtouCebxZjp5YnFs0qtF0Cgrq/kn/8O1Thxqy/Aos5lSqVGLcnaUVo3KItmBai9ZLp79w8GH/cPTzo+n0xtjZQ+4+qpx5mOw/UHMQbh7QF4K0UrPtl2TYMfnxSwvKy1JWu4zo0wJx7gv9us36tybI4wpL5953SJei66at3BKSAqO4Y2SngZ/UcsHp+jm+VPRS1vVI0kxBrHvGuMOrnc+mwI2TJ5r2+1/qSt2ah6/0/Bd3H6gxnpklRtkRiUBcH4zaRCpcHedlZOEtZqXvw+6W+1Tl1+9SiB070XY7kBFJK4LUpf6VklqK3dgJcLTcbM9bswOy1u+DxSN/4v18YmO8TbU+wVBFCIC8nCxUuD658ZzY+mLkOm/f6v2BfMZmlUr/szcokm6lwe3yzWLPvHoZ3rjwagLfsf6t6hb5mzvF68LSueOKc7gFfoGqBkki+WFMhXC7aTUOTt0VUX51duHEP7vzKn8to1fvqts8X4OxX7Mt7MJNnWGEep/SnNFYE0x3VPDAfKZL+hHpbimh1bVoLBbnZvmq5qnRYMXAKIUQuvAHdR1LKr7Srtwkhmmi3NwFg+saVUr4upSyRUpY0aBDdVm87qY2oL+1fjI27vd/HauVoJ/TXVrdfRlvROhaHKtxBq5Kp2B1+w0d/RpSLpn7uuzQtSmjxNHVSSt96t+9IJbab7FBR3zv6kNTj0d1fL8aGXYeDenjGqlVd/3dgXk4W6tfIx+Nndw+4TyrzqtRjnpneLWvj8XO6B02q7j1ciYvfnBWQAqIyVk4Ol19tRv3sW1VpNqP/++nvgx8Wbgl1d/y2sgxvTk1unumo56fhlBem4Yu5GzHbYqW+0uKLrDDJ52kM6uJQPS8HHRr5E5nHL96KVdsPBG2PXLbFu61kj7KdsWZBLto3rOEr0xpOKiv+6EHI1n1HcNEbs3Dea7/j0fHLcNQDP6PzvT/i1k8XBNz/VqXUfbrJz87Cq1PW+LaGqF9GZhUs1+/yzn5u2HUIuw5GtrWnfxu1h15BQJCekyUsP9zRuvyYYpxb0iLg8RPVfyeRVls0ptddfkxx0p5bLyVuLB5g1Uh9opIHUmwSpNhB3260fV85Pp+zATd85M8X6WOx4mhkVulSDxb17ZaDYqxGF0qkbT4ynfCeVb0FYJmU8mnlpu8AXK5dvhzAt6keWzKpk5JWuTiz1+4K+x0CBO+CSSU1r8eqCEWyTVAKZdmdZ2TcvSOE8BV/i5faL1TP6+354M94xmTFRg0u9SHph/zW9av7toAmqgm0eizWA5/zDHUR9oVppbB5z2FsTFBuqP76mE24AcBXNwzEeSWB4ysePRb/CtPyx5i+cEavZujZojam3TEk4rFlK+der5sU9zmzVzN8ef0Ay9/PiyB379xXZ+Cyt2en7Lzots8XWN7mVhYF+im7tTomMN/TDIO6OByudKOaxSrb//7Wz3fZqoTrF38/BuNuDi7Rb2w0DaT2wGG2NU7vAWSclSwdMwp1TMoup4v9hr3O4Sr0/b5mJ1Zt24/Bj0/2FYO5wPAlbfTvE6330q8pO4ixC7eEXcGIpu+L2SxjOuU1HtehQciAIZlbG60+Jz1b1DL9Nzi3j//ftjSKxrLJ9Plcb1L+sU9Mxn++CNwCXbd6Hjo1Nj8oqNXx9FW4ZkqFQb3lQw9tVS/eXpd52cH/jk6qTGazgQAuBTBUCDFf++9kAGMAnCCEWAVguPZzlbFgwx7f5exsYTlJMfzpKWEf664wqxLJpOaPFhiKZdRN4PFwRNfIUioj3VWSLGt3BJfuNzarjnVlMVepFqOfKFvNk6rHFqntXBYQaFG3Gnopzchf/jXyBvehXKj07rTaobLvcCWKR49Fx3vG412TtgjHjJmEQY9F3v8uFD2479Y0fCVZNbiItHaArm71PHxz48CoehOHO958PW8T+rSqi+cu6Gl6eySpR3qPZcD+XSPqar466dHR4vidKAzq4nC40m15gppt8g1mTEqtVZhr2tvl+Qt74VFDa4MKlwd9W9cNWBVKpf1HggNTtRm6U/z3h6Uhb8/NzsIJz/wWcJ1Vo9SvbjgG7155NI6OoF9YqPYQgPkB4aEzrAtcGH1yrfUMV6qd2bsZPrymH3q1rI0HTu2Cd644OuD2giirZUXDKldv4rLtuP7DP9H7oQlYvd1/QibNU5ZsNVTrAWRm+/5yyxUOdYKpQ8OaeOfKozHptuN81+lvMf0AE+80gF0tH6oCKeU0KaWQUnaXUvbU/hsnpdwppRwmpWwvpRwupUzPhOUEyMkSUVW6NFqw0b6VOvWz9m/DbP3YmwcZ7x6zSNuWJGOlzqwvmxWzgjXGicZrj20b0zhyA1bqPCG3M6oTW1JZqcsWIqB6rzq5EI8mSv0E9dDz0TX+SX393Knc5cED34c+/4iXXugtJ4JJ3nDvmVCrZrGIdMvmKRYN0dV+qpEEbOG2oiaasaK3GsjNWbfb9PpkYFAXh3KToG5gO29JdTV3S6fmDYSSl5MVVJGu3OVBZYoKpZgxC0oGJmH7VqpY9VQzC/oK87JNm2wW5mXj+I7WJ+CBzxf69pkmyfZtItyaCwC1YugXE6vzSpoDAE4+yrwvn14F9OsbBuKKga19jUoB4PySFkmv3mrlxyVbsetgBf43y1+syO5tS2ZCbWdesGGPZWsB9do/1u3CkI4NA0puP3tBTwxoUy+oIm+s2Iyc4pGTlYU2Sh+2qqBni9poUitxu2rUz9iq/zvJ8n7hJg1j8aASgIT7ntT7lT1/YS/fdcbtcrHm8OZkZ/n6rx6qcIcsE6+ms+hjzs3Owv4jLny3wLz3Z7TUbaVZWcJX1bFmgf/vU8+N9sZYICYW+k6JSM4T/zJZXVX1aVUXSx4cgZl3DkvI2CKdBLS63wFlYSGSPNYv5pq3oUiWiUsDVzvVz0x1ZfEm2QuIDOrioG+/HN7Zv0WifUPv0qr6xW48CYukh1mN/Bw8qDT1/mLuRsxbvydlXxBqrqCVqwa1TsFIkiNcTzVVQW42pt0xBOeVNMcT5/iToCM5SN0zytuKoDJEVGdso6Cz2kM+667gL9migtStmv7fmUdh8m3H4+WL+wSU/tcb7pp9KV9/vHeW1mxlOtXenr7Wly/58xJ/Xkq0TXGTJVzfnyWb95le/7CysmvWMuKYtvXx8bX9cWyH+ijMy8bVcX5+1Qmtn289Nq7HosyTkyVQM8aiAXZPxljluJ8WYXuQSKlb1nKzszD77mFBlQiB5G+/HP3lInwyez2u/3Cu6e366oN68prICWi9qmW4E3V1p4b+muRlC19lVV0k5zdWercM3DJ87yldMOuuYUHnA/oKozF3f9kW8+/vRNBzVhP12lfPz0HjBBV6C7dSZ5aycVI3/8SxeloRyWpXLMVc4lGjIPC7TP1MHq3k1Bnvl2gM6uKw/4gL+blZAf0o9A9yX2VLXrmh9G6k/Uwu7d8Kd53s/QLXE3sXpmjLyUsX9Q57H+OXm1NE+6Wan5OFnOwsPH5OD5yrJBlHMiNbLYIG6AeUvL/RJ/kP2FZfzI2KCnBun+YB1zU0WUlMltzsLN+MaEFuNmbfNQyr/u8kPH9hL0z693GmDVkb1vSuDqV6n7tVRdC1O7yFGBopr1ui+xvFKitLBLRgOK+kOZ45v0eI3/AaoeTiWhWGAbzFfJb+dyR6KDkmsSjR8kDP7dM8oc1+qWrxeCQ27ArOV83KEgHffaFs2Xs44LvjULm9FZetgrreERYyipRxNbxhzQK0aRD8PZXsIPfLPzdi9FeLMN4iJ12/Xg1CjRUT45Gj5dU998uqMPf0q/RtRQweR7jCcy9OWoUx45djm9Jbb969J2CuScuCnOysgOOI7rgO5rt4Tnpuasjnjke5dp6h92uN1QvKimuqqJMV397o7cN6q9KOQt3xE0lQl6w6A9/O32R6faFhwlpdzFHP5RKZc2uGQV0ctuw9gm/nb8bNw/z9mvR/PPXLbd6G3QGVFNs1jGyWKCtLxLwPPV7psKKSCFNvHxI0s2lMkq5fIz9koB1PYQ89wAl1EqIeYP5+nP/fO9RBMZoE5WRrWFSA3OwsFORmh91OlYpy5dce6+/7uPD+E0Pc01viOR0tVKqz1SzIDQjYIhFJpbB4FeRmo3TMKDxxbviAkzLXa7/9hcGPT8aqbcHN6dV85P5t6gYUCZNSQkqJv8oOYMCjk/DqFH/FvFA7H1LBrHri4gdHoGecEyVGDWoGT86Y9dw02345Y/UOHP/E5JATirEwO6F+RTumzlNyh+ItxKSKJEfMqMJtvWplFoSpnvx5JV6dsiZgQr5O9Tzf1vWiCFZb/n5c+P7DKrOJj2j102ouRLJi/H9nWufsR5rLGY3jO1rvUHvy3B7o1sz/nD1a1EbpmFEBk4XqRKVVCoLqYBJabUkpccsn801vM37O9B1Y2VkCd53cGeNvGYyP/9Y/4WMyYlAXo2cmeMvpntilEarn+0/6zWaAfl1RhnU7/fuX7WxsHKlEzrLZqUXdQlx/fFusVvIRvlf21tfMz8Gce4bjumOtv4DDbYcLRa+2+XqIvin6weexs73FcfRtHHkmZel16s6CC/uGrs6ZDvS+acneegAgoEdSqBOLSrfH17/uh5sSV9wgEdTvkbycLBTmRfa6jb15EP4xpF3IlbpkqV8jH1cNdO6WbEqOmVrT7PXKSWvpmFEAgDrV/VvW/jOiI1rV909WXfD6THS890dfgbEPZ67z3WY1U5+qnQBmVZSrJ2Ei9IQuwZM5Zjmx6vfFjNU7sHzrPtz/3RKU7jyEdRFU9a10eyzzzM3ua+WAYfJSnaSMRyT5WFcYWuXobwWzX410RdVlMXkw++7hAakHZkKtJhWPHhvUkmP3ochaKIXSsm511MjPQa8IdlGZbdHXxXPOY6V/m3pY++jJpred06e5ZZEzXfX8HNx/qrdpeqTv1USzSn0AvK2c1M/Gtv3eKuXvX9UXretXR+cmRRjQtl7Sx1g1ztxTTErp2wbQom4hjiizOW9O81eMevEi7xJ2y7qFKd/fG68GNb2rV1Yf/PYRrjamC6viHC205qGhincYg/DXLu2Dn/4ZWQ6RHpht2RNcOEfnskhuNisZr1PfTWvKQic8p4OzejXD6JM64eah7cPfOU6RrG673BLHP/ErNmn/Lt2ahS8BnUrq1hHjBMviB0dY/l7XprVwW4TbuxNtzj3DcZ920CXS6e9ls62WapXo3OysgOPkrLW7UOHy4LM53jyqTcp3qH7y9PjZ3X19Vb3Xp+Zkb4XJqmO4k9JYWBUjGtAm8OTwsLIqcdGbszDy2am+/Kpw87NSSrS/ezz++8NSSCnxv1nrsTNEc25jT7VQwcsdWoGTM3rGl2uYk2X9R8y+axiuO7ZNQOoCELq67ysRtjSw+tsKcrPDTs5v3G19zAeA6at3RDSGaJS73BEHZLULrbcBNq2dnBZa8X5G9C2OZit1K00+kzd8ZJ4DGqtDIVb/Hvx+Kdrf7e/Tes/X3uqbU1aWJXQM4TCoi8H7v/tnDNfuOIjaFqX99dKs63cdwnythG4kuWqhmCVJJ4MQAjcOaWfZR+hxpWCIU5g1el8aQX6d8aR6RNfGEfca+adWydCqTC8QWKUL8M82h9pCV6z8LbPXpn/V85zsLPz9uLa2bOs1ew+/P3NdwEliumlZ179iYXwfWOUJEqWjDbu8nzOzbUvqCkxOVlbIdh4qfSIsJ1tgkNIWwbiKVOHyBDTpdhqroO4ipT8aYH6yqa+MhptQ1neKvDujFJOWb8ddXy9Cn4cnWt7/oR+W4ZQXpqL93eOwbMu+gFXC1oZ8PyEE5t93QtxbtENV2m1YVIA7T+4cFGTdMKQt2tSvjiERVKgud7lx0GTSQW9v9HAU7YV04bbAPzJuOT74vVQZQ/xbij+atT6oKEzI+yutF1SpalfTo3kt3+pbJLK14N4s2P7NJHgatyh0X+Jk0it0prr1D4O6GKgzjiu27kfbCMoy3//dEgDx/wPrVQRTxWq84RKN09E/hrQLuk7NpVP7qU1QqvnFkxugn4BbbeMA1DLEgc8TKtH3lO5N8G8tibgji1QEGdy+vq967PtX9Q26ffW2A6keUlTUXWT6roDXLu1j+rcQpasjlW7TVS2dceZencwIRQ9Yyl2egGOUvtVT9/SElfjb+3MwY03iV0VSwer4e2LXRujfpi7eurwEAHCowntOYrb9NGxQpxzLjT33zHikxOJN+1DpljjpuakB+Y1nmfR0rV2YF3c1xuwYcuraNayJSbcdjzoRFKY46+UZ6Hr/T5a3q1v6I6UXkgrl3m+X+C4nOvcxEnYXu/vk2gG4Moot+3pwH0lOnW7vocRVjI/0bSyl9J3XndzNeptrMjCoi0GxUiXPONuvl7DXGbcvxlOV6I3LSmL+3VipxwN1z3qiytymUvM6wVsK1FVWtbpSpMVswtG3QhypNA/qKt0eLNq0B4B/pU5fUQw10yeEwAit3K87xRUlneCDq/vhcu39Wt1kZWvFtv0RtRaxS3ul5LZ+0jWia2Mcm8ZjJjJ60KTZstUc2YFyV8S9zPTjUruGNQICH+NX4Y+LtwAA9iWgFdCasgPacwQ+Sb/WdTHu5sFxP74ZPehtVS8w2M3PycYn1w7AMW29q5Rz1+2GxyNNq2Of+uI0zFUKmBhFO0E7dVVggKy+tmbftYmQneT0lVC5UkBsLQLMetuGsjkBO0eOLq6DY6LI20pWhchIRbtzR/+su00myXtpRc+MlTv//fmChFWHfeC7yJrHv/abv4ZCo1qpzW9nUBcDdc9yM8PeY+OsQ+cmgasoZr04wtF7qpzQpVGYeyaeuu9a7WsV6YxqOjELfmoa+ru9fUUJXrqoN4QQ+PL6Abj3lPhyhPK1LSFWjU8fHbccd3zp3Xut5/W9cVkJXrm4d9DYjPQvOLuShp0u3rLPyfSksl0p3nwUIjv8tGQrPp69Puh6q5WjklZ1UJCbjS4hKu/puWMeXxEMgYFt/cfUa96fA8Db8Pna9+egVCsSEioXJhJfz9uIYU9NwbMTV2LV9sBV/sfO7o4uTRNfLVD32XUD8OX1x5jepn+HfTN/M056bioe+H5J0H32H3Hh7Fdm+AJcoz1KUBZLgbRBj02O+neiZXbszs4SprtvzNx9cufwd9KYBbmx7LCKNn/s2YmRt2uw4pHR9WcLVUsgHel/25LN+4Jy6G7/YiEAf09D3cRl23Dai9MT8vxm1bJfvSQ4per3Nf4dA7kh8kGTwVn/omnihUneD99J3Rpj0m3HBdxm/PDrJe31vnWxfIh++uexllWDku3U7v6VRrWscDKSwpOtpFXdoOuGdw7cbz+0UyOM0v7mPq3qxt+gWZsAsJopnV2qfPi117dhUQFOClGZShcqz4ACjTnrqKDrUtXzMRZqQG+V10qUzq77wLxIwfhbzFe19G3u+45Yr6o9oKUxeLSTfCGAk45q4itKBngnuf43az1+VnLpauTnYPnWfVi8aW9EPa6M9KbXz05cFVSAI5Zy+9Ho27quZTVb9Ti8Ytt+zFu/x/JxXpi0Gocr3AHVnwHgnq8X+y6fqpTCz9PyoG9IccqHGbdJAZw1j5wccVGov4Wobm102GQCINYCPN9o/dYicbEhTzIWHinhlFOzBfeFbjdk5rM5GwB483NP1PIddXrBOLP2Holq+N6xceDkzaijmmBop8DFliOV7oDiKLkpaC+kYlAXgwXayeANx7fzBW0zRg/F1NuHBN1XP1GfXRp7MQshhG1BlPq8qU74TDSz7YyRloqPlRrEm1VBUr9/ot3i0aJOIa44phivXdon5vFlCrPAaMveIyb3TD8X9Wtl9xCIEsZYWe/9q/oG9MwKVVXwU+2kDspKHQDs2O+v1vjBzHVBWzzfmPoXRj47Fae8MA2P/bg86jFPX+2ffNtxIP7S83bIEgJXv/cHbvp4Hp5Xmnj/ruQhvqVU765we5CXk4XbIyzOlszemMYcqudjaI59du/mOK+kecB1Ho8M2E7r8UgcrAgumBLr39ajub+q8rx7Twh53yd/XhnTc6iiXakD/EVgkt0UW/fh1f3w1uUlqGVRYDCU848O375p057DSSsmptZG+PGfg/HSxb2D3hsXvzkr4OdUT74zqIuDmlfWtHY1X3l8ldOWt0MpKsjFNzcOxCNnBq96UHhmzXfVA0q0H/6sLIEHTuuK9iyUEpY+e9fJpGrpbSd2SPVwouL0yRQilXHy6tgODXCxMnHRv03wjgojfaVO/2iojYbv/25J0GrFH6X+nRJm1TAXb9qL4tFjTb+jgcDvjYvfnBlwm1NSmt0eiRnatrD/zQreFmsmL4pVyCeSWBFbXV09pm29iJprGz11Xg88fk4PX7782IVb0Oaucfj3Z/7iMPuPuEyrUObG+B2sToqnYsVGSmmZs2pFb091urbNP9pcwGgNal8fwzrHlkoUSVHCbCGC2v4kKnfQpazYWq34G3dlxVskKFqsjR2ly9+e7bscyexNVdoiV5CbjZ4taqNni9p2D8WRzE7O1Ya8oXrxUHwa1PRuXzqrdzM8Mi5wpv4fKeidF4ufbz0WC7RWKERVRbgTrEv7F+PDmaGDjj1aRTt9VcK4OhFqtcLsZGzcIm++2XitqbhxoqxDo5pYvtUb8Om/3rVpEZZs3mfZ0ijdLNvq34J2mpKnm5eTZVksJZpcxGSevKo7LeKd5Fqt5UTqTe2/mrfJd1uP//5s+jsFCWjFk4pzwYUb96JF3eh6zPVrUw/jbh6Mzk1qonX96hjZLbjhfbo4bFEhVK0r0MSkx16lW+LWT+fjgVO7xrRCqFNX6tSCiaGwpUGaU7fQFRWEj4mrwiz7M+f3wF0np6Y/XlVmdtBTD5pWX1gUv4Y1C7D4wRH42+A2aV0gRdWhUU2cWxJ+uwmRk4RLJejYuCaGh5jJr3B5fKX3t+/3bqGO5jPtMVla03fUPD1hJU545regFTuznpb/GNIOfz1yctiCVulC/bPVrXaX9bfe3v1yhE26geQGLY1rFaBEC+yM1TdjtWhTZDnVA9vVQ68ETGQX5Gb7WlBYSUSrKL0vZDS6NC2CEAKXDShGw5rpW9nc6i32ntLv71itb6VxV87X8zbh5Smr43p+fRvw8odGJq3Sa7yccXaThk7v2TSiPLdoGkGmqzN7Nce1x9qfLJ1IF/ZtgeopboStzwZbibScN8WmRn4OhBB44pz4GuESUWjxnpy+eXkJ3r7CfwLcsKa/UMiDSoVHvYDFhX0Di0yEOjabVQs2BiR6+wJdB5Mt7jnZWXH1ME0Fq2BXDWwTVTw52dvMkpmzF8rlA4oTVtMg3LbDeErvz14be90Gp6hdzTzvT61Iq/9bnWnSM/FInFVw9e2X6bwDj0FdFPYrVbkiLbCxRZnhe/b8nokeEsXo0bO6Y8l/R6bkufTKpTPW7MS89btxzzeLfI1G1apmHU3yvSjxqsLqOVE6O+KKf9eBWlXu23/4qwh+pOSD6U2hQxVXicTWfYFFk56ZEFhe3myLZTqf2Oms+qM+/uMK32Wznl86s3LtVpJdBXTltgPh75QEyQ5WFz7grwIZT1B33mu/AwC6K8VZqpqW9cxbac0t3R10ndp3WGc2gTHrr52495vFwTeY0FfqjOcQVjmekeQHJxqDuij0eXii73JhhKs86uzSGSYzB1T1qbN8Z748Ax/OXI+v523CzL92YseB8hC/ScmwfZ8zql4SOdWRBG8lN+Yb3zTU25+sd0vzdh9zDNWm1X6ym02q3hoLh6wwbL8069+WLpND1eIMaCtDLNUZq06GkuyccLuOlcn+d1bzPxOx/TKdW/Uki9kkkl4ARvXBzHVB153/+kx8MHNdQNE6K26PBzlZwdXon7+wl2mhoA+v7hf2MRONQV0U1A9cpEFdrP1NqGor3XkQF7w+M/wdKeHUg3Q0TWmJKDLlJitE953SBbPuGhbT4zWoGdinTS+3b7UtTi94ojPLiVNdNsA8r2z19gMoHj0Wr0wJzi2LtnR8skwfPTSu33eFWB0a0TWwaIbZibIuURUG002iV2RLx4yyvG1BBgZk8dC3Uq/beSjoNiFEyNfa6A+T1T4jl1taBvlm19tR/Z5BXYwinZW69tg2aFAzn73E0sR9p3RBr5a17R4GVhu2kqTLrG8m6KXM7kfTlJaIIlNuMnN+1aDWaBRHufSzevt3urw25a+ErGrozNIpRjzzG6au8hZGM3uueLbKJVLd6nkY0rGB6W2DtaIRZjZolZddFhPPg9rVD9p6ePOwwErB6mHLKe2bXo/yXEyvehqrNy4rwQOndrG8Xc1v/Nv7c+J6LgBoXie66pdOo+ba6q2KrN7/0dC3r4ZS4faYrtoDwNHFgVstf/tPcN/qVHDGpzBNqP2sIv3+alq7Gv64e3jQjBfZ46pBrfH1DQPD3zHB9Aaful+Wbw/4efEDgX1VKHm6NauFu07uhMm3HW/3UIiqJKtcrmhdd2wbdG1a5P0hzKaXufcMj+qxb/p4Hr6et1F76OAHX7FtP9SpNmPRkWi2JibbyxebByp3jOxkuVrx9w/nAgguCqMzqxJar0ZgoQp1tTLZOYZHNfPmio06qklcj9O/bb3onjfOHLUTujTCFQNbB1w37ubBvsbkBTmJKdimvz8/ubZ/Qh4vXam5tr7+s02KkJedFXZlLtp2D0aV2g4B88cOzPezyv9LNgZ1UVCbUu4vd9k4EnKaUDOmAFAtxZU4M921x7ZF6/qR9ZkhougYc+oW3HeixT1Du/Pkzhh782AAYWM61KuRj4HtIj9h/37BZtz6qbc1wpdzN4W5d/DunF9XbLe4Z+oZjx/z7zsBpWNGoZsWCJntTjmgncNYbfkzKw5Tu1oext48CEO1IhRqUJfsgiKPnnUUTj6qMZ46L77qxUVRtqBQ2z8kSpemRfj2H4NQOmZUQqp6bth1yDeR0ryOPcFEKt15krfF1jvTSjFlZRmWb9nnC/BCCbfDbsve0Nu0K1yelDcTj1Z6jy7NqFswZq7ZaeNIyGlCfXGbld4lInIqdaXu3lO6xNXwVxdJIYPRIwNzZP9+nHkrnj/X+/NnHvtxuWURjge+X+q7bNwiP9Skul46WPrfEahdGBiIzFu/J+h+ZnlIugdO7YJHz/QXfvjxn4Nxao+m6NS4Jro2rYWXLvJWxWyr5NglO4OgW7NaePniPnFXOo2W1Xa7dHLai9PsHkJK6a3Cnpm4Epe/PRuTV5RZ3re+srpcM0xv6aFPTgl5e4XLeqUuXaT36NKMulJ3JwssUBTyQ2yxeIatLoioCtFX6r65cSCuGlickMc8FEGPKeNJ2+iTOuH0nsHlxq9+9w/f5VcibLCdmy3QrVmR7+dwPcfsEmm7pVCuGNg6IBDv1LgIL1zYy9eXr1peNt64rAQfXN3Xd5/02YwavasM2yNVqViZGR7ne2n3ocrwd6pCFmzYE/F9HznzKN/lcL2AD4ep2lvpliGDuhlxFi1KBAZ1USh3eVCveh6+uuEY9G8T3b5symxWXwRmZXCJiJxMLzFeLTc7YY2bf166Lex9zHqlmZ2Ux3ISnJ0l8MNNg1GnMBf1krAlL15/G9zan38YwjDDCmPzOtUCitBE6oQujQL6rNYOc8KcTvTdMbWq5WLJgyNwn0khE32CwCy3MNGGd07PVd90tX6X9Sqz0cB2/tQXq6JAkSoPs/2yae1qWHD/ibYGdwzqolDucqMgN9uyNw6RFastHOf0aZ7ikRARJZe+/dJYXCQes+8O3w7BLGfGeBIW65i27fNu0fzj7uExt2ZIprtHdfHlH4bST2mIvONAOTbuPpyQdJKGcVQ2TbWnz+uB5Q+NxKy7hqF6vjd4++Ta/ujSpAhXD2qNxQ+OwMNndENeTlZScuqMzj+6BQCgY6OaAdePGb8cxaPHBlVanbKyDK9NWWNZ5KaqUycTwslXJtTdcRY3ClUoRVerWi6a1ravAimDuihMW7UjbL8bIjNmPXzeu6pvwmaxiYjShd7SIJH5Tw1rFuDJc0MXyTD7OtW/e1vWLcRVA1sjN84m2TnZWY4p36+77jh/65ZTe/i3oz7w3RIA5g3ZqzIhBApyswPen/3b1MO4Wwbj3lO6oEZ+Dk7v2QwrHz4pJTl8+nmAsen9q1p/xE//2AAA2Lj7EIpHj8Xlb8/Go+OXY9hToXPAqqp7T7FuEWGkflYrPYHBcSR5urq/yg5gysoyrN4WX4uLZHPWN5PNtu83T6YmCscseIu0gT0RkZP4VuoSVK5dt+dQRcjb1ZwZfSviuEVbAHi3bOXlZAXkxmeKS/v7m6urq5k/LNwS92OP6t4ELetW/YqLqbbroP+9/tOSrQCAz7TgzkhfPTL2Eayq+rauG/5OJg5XuOFSVj0/sXg9zQzVAuiDEeT22olBXZTYJJoSJdl9fYiI7KAXSslP4PZLwF+G30pBbja++0dgH9IdB/wnx3k5WahwezCoXegWM1WNWuY+J0sk9O9/6aLe+O12exotV0X66lHvhyb4rpu6agcA69z8s3t70zjOK2E6RyjLt+7HKS/4K4Xe+dUiG0eTHAzqohTvnlwiHbdeElFVVF7phhCB+SyJoG6Fm25RjEBoLcPNvl4/n+OdmS/deTCh43KSnGwRtNLx9hUlNo2GjNbusH5vWp0zNCry5pg1qWVfLpdTLN+6Hzd9PM/uYSQNgzqiFGtay5tQ3k7p8UNE5ETfLdiMzYZc8+cnrYaUiZ+4Unc3NLMoRtCqvndV6vrj2gXdtkXLHdu4O/rc+B7Na0X9O+koNzsraMWnkVLkRA8QyB5Dn5riW+k2mm9Syv/iN2fi2YmrAGT2TrI3L4t8YuL7BZsD+k4fXewvfhhNnl06YlAXIY+2QndLhuxZpuSZcecwlI4ZhRr58fcTIiKyi9sjcfPH83Duq7+n5PlWbA1fpKCoIBelY0ZhVPcmUT12uB5Wp3QP7nfnRDlZIqgas5pn16cVq3vb4QGlrcItn5ivJE0waesxfXX8lUurguFdouv1d/GbM32XL1FyTp1eDJFBXYT0RNVUVEKiqunHfw4OyvcgInIqPR1h2z5/9URPElMUYily8vAZ3XyX21vsjhjeuSFmjB6KYZ0a4vc7zbd1FlWrGpNw2VkiaKWuuL4/5+70ntH3rKP45SvnlpOXl9k4kqpHbUCu+6N0t+9ypdK/zulpMQzqInT9R38CSGzfHcosnRoXoXvz2nYPg4goIcxyzP/+4dykPd8FWj+vaKhVhod1Np/N79ykCNXzc/DWFUejSa1quO3EDgCA20d2xNTbh2Bw+/o4o1fVCHaECF6py1eqlI7o2jjVQyIE9rI1a4E0flH8lUoz1UX9Woa8vcLlQRMtLWbxpr0h7/vqJb0TNq5kqBpTTynElToiIiLApfV9cinB3c8mW8QS5ZgYqjaqE2mjjmri6/2lMqbR/GNoe1wxsLVvi/wHV/eL+nnTmVkVxRmjh5oWl6HU8ChvQrPcuNu/XJjK4TjGnSd1Ctg+aeXcPs3x+dyNprdVuNy+lmXXfTAXpWNGWT7OyG7RbetONS47RWnj7kN2D4GIiBxOCPG2EGK7EGKxcl0PIcTvQohFQojvhRBFdo4xHI8DWr61a1gDp/ZoinP6NLcMWiSCVxyrWs5zd6XQi1lQ17R2NVZPtNF3Czb7LpvV6th/xN/OQ19VIuDKga1RPYLP6nXHtbW8rUZBruMLpOhsCeqEELcKIZYIIRYLIT4WQqT1O1TdYrL3cKWNIyEioiriXQAjDde9CWC0lPIoAF8D+E+qBxUNlw1R3bBODfHcBT2j+p0XLuyFJ8/tEbAaompqUUmzKvnq+mOw4mHv201t6XD5gPCrHJR8RUqhnnC5o69c0ofVszVmW1XNhHq9zurVDCeEKLTipIAv5UGdEKIZgJsBlEgpuwHIBnBBqscRjeOfnOy7rFaJIiIiioWU8jcAuwxXdwDwm3Z5AoCzUzqoKK3YZl2NUq3ml0hvXXF0zMU8jCd2U28fglcv6Y2L+obOuakKcrKzfLlzajn303pWjaqeTnfj8f4WHBXu0EFdq7qFuGdU52QPyRESUdgkK0uErPqqL+zc7IDq93ZFKDkAqgkhcgAUAtgc5v622rDLX+KUX4BERJQkSwCcrl0+F4BpZRAhxLVCiDlCiDllZfZVyrv54/kBP+84UO67fNmA4tQOJgKFeTl4/sJeAIBOjWuiRd1CjOzWxPEV76KVqxTlyOZEdVro0Ch4Jal+jTzT+woBkw3DFAt9m/VF/cxXrNeUHUC7u8cDcEahxJSPUEq5CcCTANYD2AJgr5TyZ+P90uWgpercpAi9W7KHCxERJcVVAG4QQswFUBNAhdmdpJSvSylLpJQlDRo0SOkAVflKbpbbIzFaKeaQlaaNkKtr1TAzYcullbYN/AFETpr+O2WanOzg0/F7T+mCphb5c0cX1032kNJa6/rV436MoZ0aYtzNgwFY59AOe2qK73KuAyZA7Nh+WQfemcjWAJoCqC6EuMR4v3Q5aKmObR995S0iIqJISCmXSylPlFL2AfAxgOBSjWlEPRGqdHswZ93uEPdOD8e0rY+B7erhrpM72T0U2wzp5D+nMmtLQekhPycbT53XM+j62oV5qJGfg7YN4g9snOqbGwZi4r+Ojep3vvj7APxdKZhyzaDWaFmv0PL+L/+6OuBns6qk6caOsHM4gLVSyjIpZSWArwAcY8M4onb7yMw9CBARUXIJIRpq/88CcA+AV+0dUWhHt/bvXKl0exwRIFTLy8ZH1/RHu4Y17R6KbdQVhya107pOXUbLyxGoZ9iC+e2NA32XaxZ4i6uc1K0x3risJKVjs1utwtyoP8MlxXUx+iTlPN4iRvtzvXdy6vEfVwRcnxNhURY72VGzdz2A/kKIQgCHAQwDMMeGcUTNCVE6ERGlPyHExwCOB1BfCLERwP0AagghbtTu8hWAd2waXkRqKRX7Kt2SiT4OoW6NrVtonrdFqVeQm4Ujlf4iKbnZWUHtDXq0qO27rKeCXjO4Nfq0yuztmLGwKt67fMt+0+2YWQ7IvU15UCelnCWE+ALAnwBcAOYBeD3V44jUhl3sS0dERIklpbzQ4qbnUjqQMMpdbng83hUuI5fbf8b5zvS1AU3IyRnMcrnIHmpABwDZQli24QD8C00OqrifVowtWarnZeNghRuVbg8WbdwbdH8n5J/a8mmWUt4vpewkpewmpbxUSlke/rfsURmmtCwREVFVNeSJX9H5vh9Nb1NLr78waTVqFFStht1V2WfXDUha2wlKjHkb9oQsCKJXbWVMF50hHb05pcbt4i9d3BsA0LJeIY643EG/54TdepyiCWPXQdPiY0RERFXe5r1HLG9TV+oAoGx/2s7PkkHf1nVxxcDWdg+DQjhc4UZBrn+F/LVL+wTcfm6f5gC8fesocnprsvaGnLwGNfMBAOWVHhyuCA7qmFPncEcq3Tjn1d8BAP93ZjebR0NERJQe3v+9FB/MXGd6W98ML7dOFIvW9atj7Y6Dvp+NK0MjujYO+PmCvi1x/tEtMq7PYrzO7NUcp3RvGtCvEfBWGwW8OxBamgTKwqqyShrhSl0Im/b4m44v27LPxpEQERHZx3gMvO/bJZb3ffeqo5M9HKIq57PrBuB/1/Tz/Xxsh/BttBjQxcYY0AH+vpvllW4s3hx8zl/hgHQsBnUhPDJ2me9yYR4XNYmIKDPtOBD51koeL4mi16BmPo5pVx9LHhyB967qy4qWKaYHdRVuD57/ZRUABPQCnLx8uy3jigaDuhAqlSTKE7s0snEkRERE9jFW5iOi5Kien4PjOjQIf0dKKH375d1fL/Zd98NNg32X1zugGj6DuhDURe0uTYtsGwcREVGqqatz5SbV4IiIqor83OCQqFpeNo7XqmUe7YBcYQZ1IUxZWea7rEfwREREmeCRcf4UBGOlSysOqPpNRBQkz6JnY0mrOgCARkUFqRxOTLjxPUJO6E9BRESUKBUu/5bLhkX5vsvPTlxp+TvsP06UWBNuPZa96FIgy+I8/6pBrbHrYCUuG9AqxSOKHlfqTMzfsAfFo8faPQwiIiLbnNStie+yWvzk2YmrLH+HE6BEidW+UU10aFQz/B0pKQrzcnDfqV1QPT/918EY1Jk446Xpdg+BiIjIVrlKs123x7xQSs8WtQN+dnOpjojIFgzqIlA6ZpTdQyAiIkoptQeWVYumxkUFGHPWUSkaERERWUn/tUQiIiJKOXV1zmWxUvfYOd1xqMLl+/ms3s2SPi4iIgrGlToTg9rVt3sIREREtvp+4Rbf5U//2GB6n1rVciGUBkBPn9cz2cMiIkqKutXz7B5CXBjUmVATvS/p39LGkRAREdljrBLUfTt/MwDAw5w5Iqqi5tw9HN/cONDuYcSMQZ0JNdF7ULsGNo6EiIgofbz62xrf5RO6NAIANK6V/v2biIjCycoS2Hu40u5hxIxBnYlyl9t3eWS3xjaOhIiIKPWsVuTmlO72XVarYxIRVQVO7srCoM7EH8pBi4iIKNNUWhRGyVHOeCrd3IpJRFVLt6a17B5CzFj9koiIiAJY9Zubv2GP7/KEpdt8l/u2rouigtxkD4uIKKlqVctF7cJcjB7Zye6hRI1BnYF6IPvf3/rZOBIiIiJ7mK3C7TtSie37y03v/9l1A5I9JCKipMvKEph/34l2DyMm3H5poJZtPqYtWxsQEVHmcZl0G7/hwz9tGAkREUWCQZ2BW3pnJwe2q2fzSIiIiOxR7vIGdQ+d0c133bTVO+waDhERhcGgzqBBDW/jwTtP6mzzSIiIiOzx2hRv64KxCzdb3qdjo5qpGg4REYXBoM7gYLm3nUFhXrbNIyEiIrLH+l2HAPiPiWZuGtYuVcMhIqIwGNQZ/PvzBQCAglwGdURElJn0YK5Gfg7uPtm7c6Vr06KA+1gUyCQiIhswqLNQjUEdERFlqNmluwAAHRvXxN+ObYPcbIGmtasBAG44vi36tKqDIR0b2DlEIiJSsKWBhdqF7LdDRESZ7cYh3i2WlW7p60vXrmEN3O7AHk5ERFUZV+oMsrMEAEAIYfNIiIiI0o9+nCQiovTBlTqD9g1roEXdQruHQUREZDuJ4MQ5BnVEROmHK3UGFW4P8nP4shARETWokR90XQ6DOiKitMPoxaDC5UEegzoiIiLTVITsLB4jiYjSDbdfKqSU2Lj7MBoVHbJ7KERERLayKhiWzZiOiCjt8KtZUba/HAAwd91um0dCRERkn9qFuTi9R1Pfz6d0b+K7zJU6IqL0w29mxaY9h+0eAhERka2klNhzqBI7Dlb4rvtD61sHANmsDk1ElHYY1CnOfHmG3UMgIqIMIIR4WwixXQixWLmupxBiphBivhBijhCirx1jW7xpHwBg7MItvuu27Sv3XT5S6U75mIiIKLSIgzohRCshxHDtcjUhRM3kDYuIiKhKexfASMN1jwN4UErZE8B92s8pZ9bGoFuzIt/l7Gyu1BERpZuIgjohxN8AfAHgNe2q5gC+SdKYbPevEzrYPQQiIqrCpJS/AdhlvBqAHj3VArA5pYPS3PHloqDr9NU7ADCJ+YiIyGaRrtTdCGAggH0AIKVcBaBhsgZll+7NawEAbhrazuaREBFRBvongCeEEBsAPAngTrM7CSGu1bZnzikrK0v4IJZt2Rfydm6/JCJKP5EGdeVSSl/GtBAiB1Vwrm7hxr0AzPvyEBERJdn1AG6VUrYAcCuAt8zuJKV8XUpZIqUsadCgQdIGo265VPVtXTdpz0lERLGJNKibIoS4C0A1IcQJAD4H8H3yhkVERJRxLgfwlXb5cwC2FErRvXBhb9/ls3o1812uVyPfjuEQEVEIkQZ1owGUAVgE4DoA4wDck6xB2WHvoUq7h0BERJltM4DjtMtDAayycSxoXb+67/Kwzo0AAP8Z0dGu4RARUQg5kdxJSukB8Ib2X5XU478/2z0EIiLKEEKIjwEcD6C+EGIjgPsB/A3Ac1qKwxEA19oxtg6NaqBtgxoB153UrTGev7AXRh3VxOK3iIjITiGDOiHEIoTInZNSdk/4iGx2QpdGdg+BiIiqOCnlhRY39UnpQExUuDzIywncyJOVJXBaj6Y2jYiIiMIJt1J3SkpGkUbqM1eAiIgyWLnLg/yciNvYEhFRGgj5rS2lXBfqv1QNMpUGtK1n9xCIiIhsIaXElr1H8NmcjXYPhYiIohBp8/H+Qog/hBAHhBAVQgi3ECJ0IxuHGd65IeoU5nJ7CRERZawjlR4AQE4WW/sQETlJRIVSALwI4AJ4SyyXALgMQIdkDcoOE5dtt3sIREREtvryT+8K3bklLWweCRERRSPiTfNSytUAsqWUbinlOwBGJm9YRERElGr3fLMYAPDx7PU2j4SIiKIR6UrdISFEHoD5QojHAWxBFAEhERERERERJUekgdml2n3/AeAggBYAzk7WoOxy+YBWdg+BiIjINp0a1wQAvHxxb5tHQkRE0Yh0pW4HgAop5READwohsgFUmdr/UkoIARRVy7V7KERERLY5sUsjLN+6HyO7NrZ7KEREFIVIV+p+AVCo/FwNwMTED8celW4JKcG+PERElNEOVbhRPS8bWax+SUTkKJFGMQVSygP6D9rlwhD3d5RylxsAkJ+TbfNIiIiI7HOwwo3C/Eg38RARUbqINKg7KITwbbAXQpQAOJycIaWWxyMxe+0uAN4ZSiIioky1puwAyvaX2z0MIiKKUqRB3S0APhdCTBVCTAXwCbxFU2IihKgthPhCCLFcCLFMCDEg1seK1wuTVuPq9+YAAD6bs8GuYRAREdlOn+QkIiJniXSPRWsAvQC0BHAWgH4AZBzP+xyAH6WU52itEmzbyvnJH/5ePA2LqkztFyIiIiIiyhCRrtTdK6XcB6A2gCEAXgbwSixPKISoBeBYAG8BgJSyQkq5J5bHSoRsJRn8nlGd7RoGERGR7bo2LUKb+tXtHgYREUUp0qBOTzYbBeANKeVYAHkxPmdrAGUA3hFCzBNCvCmECDqCCCGuFULMEULMKSsri/Gpwtu4258a2KFRzaQ9DxERUbqrlpuNJrUL7B4GERFFKdKgbpMQ4jUA5wMYJ4TIj+J3jXIA9AbwipSyF7zNzEcb7ySlfF1KWSKlLGnQoEGMTxWdarmsfklERJmr0u1BThbb+xAROU2k39znAfgJwAhtq2RdAP+J8Tk3AtgopZyl/fwFvEGeLUYd1cR3OSebBzIiIspclW6JXB4LiYgcJ6JCKVLKQwC+Un7eAmBLLE8opdwqhNgghOgopVwBYBiApbE8ViI0KvJuMxl78yC7hkBERJQWKt0e5Gaz8TgRkdPY1WH0JgAfaZUv/wJwpU3jQKXbgzqFuejatJZdQyAiIkoLLg9X6oiInMiWoE5KOR9AiR3PbVTh8vAARkREBO8xMYcrdUREjpPx0Uylh0EdERHRwXIX9hyqQB6PiUREjmPX9su04fZI5g8QEVHG6/XfCahwe+DySLuHQkREUcr46TiXW7LqJRERZbwKtwcAMG/9bptHQkRE0cr4aMbbk4crdURERAAwomtju4dARERRyvigzu2RTAonIiLS7DlcafcQiIgoShkf1FV6JHKyMv5lICIiAgCs33nI7iEQEVGUMr5Qyp5DFfBIJoUTEREBQMOifLuHQEREUcr4oG7hxr12D4GIiChtCDAlgYjIaTJ63+EzE1baPQQiIqK0kpfDoI6IyGkyOqh77pdVdg+BiIgorTAjgYjIeTI6qNNls6UBERERAG9VaCIicpaMDupqVcsFAJzUjT15iIiIAAZ1REROlNFB3V6tFw+3mhARUaarnpcNALhiYLG9AyEioqhlfPVLAKhfI8/uIRAREdmqXo18nNiqDro3r233UIiIKEoZvVKn+9eJHe0eAhERZRAhxNtCiO1CiMXKdZ8KIeZr/5UKIeanckxuj4RgijkRkSNl/Erd0cV1fLl1REREKfIugBcBvK9fIaU8X78shHgKQEobqXqkRDajOiIiR8rYoM6jJYL/Ubrb5pEQEVGmkVL+JoQoNrtNCCEAnAdgaCrH5PZIVoMmInKojN1+WeH2AAAu7NvS5pEQEREFGAxgm5Qypc1UPRIQXKkjInKkjA3qyiu9QV3bBtVtHgkREVGACwF8bHWjEOJaIcQcIcScsrKyhD2pR0pkZ+xZARGRs2Xs1/eXf24EADz3S0onQomIiCwJIXIAnAXgU6v7SClfl1KWSClLGjRokLDndnuYU0dE5FQZG9Qt3uzNPz+9Z1ObR0JEROQzHMByKeXGVD+xR0pkMaeOiMiRMjao++rPTQCAi/u1snkkRESUaYQQHwP4HUBHIcRGIcTV2k0XIMTWy2Qqd3mQl5OxpwVERI6WsdUviYiI7CKlvNDi+itSPBQAQLnLjQqXB0UFbPFDROREGTsl17RWAQCgZd1Cm0dCRERkrwNHXACAGvmc6yUicqKMDeqO79QQ9WvkozoPYERElOH2a0FdzQIeE4mInChjg7pKlwd52UwIJyIiOlDOlToiIifL3KDO7UEuE8KJiIiwcfdhAOBxkYjIoTJ2Su6b+ZvtHgIREVFa0NvTNayZb+9AiIgoJpySIyIiynCVbg8AIDebpwVERE6Ukd/eLu3gRURERAzqiIicLiO/vQ9WuO0eAhERUdqodEsAQC4LiBEROVJGBnVSSruHQERElDa4UkdE5GwZ+e2tz0jedmIHm0dCRERkv0oXgzoiIifLyG9vl8d78KpXg1W+iIiIXB5uvyQicrKMDOrK9pcDALgLk4iICKjg9ksiIkfLyG/vh39YBgAYv3iLzSMhIiKy37z1ewAAOVlcqSMicqKMDOqWb90HABCCBy8iIqLmdaoBAHK4UkdE5EgZ+e2974gLAKtgEhERAd50hJoFOXYPg4iIYpSRQZ3utB5N7R4CERGR7dweyXw6IiIHy+hv8OM6NrB7CERERLZzeTzIZj4dEZFjZXRQJ8ADGBERUaVbIpdBHRGRY2V0UFeYl233EIiIiGzncntYJIWIyMEyMiu6U+OaKKqWi+r5GfnnExERBXB5JNsZEBE5WEZOy1W4PWhYM9/uYRAREaUFl1siJ5tBHRGRU2VkUFde6UFBLrdeEhERAd5CKTlZGXlKQERUJWTkN/iRSjcKcjPyTyciIgri8kjkcqWOiMixMjKyOVzpRkEOV+qIiIgA7/ZLtjQgInKujAvqpJTaSh2DOiIiIgCoZPVLIiJHy7hv8MOVbngkcKjCbfdQiIiI0gK3XxIROVvGBXVfz9sEAHh7+lqbR0JERJQeXB6JLMGgjojIqTIuqLv/2yUAgMsHtLJ5JEREROnB45HI5fZLIiLHyrhvcJdHAgBO79XM5pEQERGlB67UERE5W8YFdTqXW9o9BCIiorTg8UjksPolEZFj2RbUCSGyhRDzhBA/2PH8TWoV2PG0REREacfl8bClARGRg9m5UncLgGV2PXmLuoV2PTUREVFacXskclj9kojIsWwJ6oQQzQGMAvBmqp+7Qc18DGxXL9VPS0RElLZcHjYfJyJyMrtW6p4FcDsAj9UdhBDXCiHmCCHmlJWVJeyJ87Kz0KiIWy+JiIh0Ljdz6oiInCzlQZ0Q4hQA26WUc0PdT0r5upSyREpZ0qBBg4Q9f4Xbg/ycjK0PQ0REFMTlkchhSwMiIsey4xt8IIDThBClAD4BMFQI8WGqnrzC5UEeD1xERGQjIcTbQojtQojFhutvEkIsF0IsEUI8nqrxuD0ertQRETlYyqMbKeWdUsrmUspiABcAmCSlvCRVz1/p9rDBKhER2e1dACPVK4QQQwCcDqCHlLIrgCdTNRjm1BEROVvGRTcVLg/yuP2SiIhsJKX8DcAuw9XXAxgjpSzX7rM9VeNhTh0RkbPZGt1IKX+VUp6SqufzeCRcHsmVOiIiSkcdAAwWQswSQkwRQhxtdqdkFBJzeySys3hsJCJyqoz6Bq9we4ttcqWOiIjSUA6AugD6A/gPgM+EEEHLZ8koJObyeJDLPnVERI6VUdFNpRbUsfolERGloY0AvpJes+Ft+1M/2U/q8Uh4JJhTR0TkYBkV3VS4vEEdt18SEVEa+gbAEAAQQnQAkAdgR7Kf1C0lADCnjojIwXLsHkAq6dsvGdQREZGdhBAfAzgeQH0hxEYA9wN4G8DbWpuDCgCXS6lFXEnkcnufgjl1RETOlVFBXXmlN6gryOWBi4iI7COlvNDippS1+NG5PN5jI1fqiIicK6Oim3KXnlOXbfNIiIiI0oPbo22/ZKEUIiLHyrCgzg2AhVKIiIh0Lg9z6oiInC6jopsdB8oBAPncfklERATAv1LHnDoiIufKyG/w7OC2P0RERBnJ5QvqbB4IERHFLKO+wvUKX0XVcm0eCRERUXpwu/Xtlxl1SkBEVKVk1Dc4k8GJiIgC+apf8thIRORYGRXUMRmciIgokD7hmcXUBCIix8qooM4jmQxORESkcktOeBIROV1GRTcuNw9cREREKv3YmM1jIxGRY2VUUOcv28wDFxEREcB8cyKiqiCjgjrm1BEREQVysU8dEZHjZdQ3uFur8MWVOiIiIi/fLhYWSiEicqyMCuoq2IuHiIgoAFMTiIicL6Oim4d+WAoAyGbeABEREQDm1BERVQUZFdTpmFNHRETk5WJqAhGR42VkUMcDFxERkdfew5UAOOFJRORkGRnU8cBFRETkdceXCwEAa8oO2DwSIiKKVUYFdUUFOahZkAPBCl9EREQAgHrV8wEAtQvzbB4JERHFKqOCuvo183FchwZ2D4OIiChtXDWoNQCgV4va9g6EiIhillFBXaXbg7zsjPqTiYiIQvKwpQERkeNlVIRT6ZLIZVBHRETk4/KwhysRkdNl1Df4EZebfXiIiIgUbrY0ICJyvIwJ6o5UurHnUCU+mrXe7qEQERGlDf9KHYM6IiKnypigbv8Rl91DICIiSjtuj0SWALIY1BEROVbGBHVSSruHQERElHZcHsl8OiIih8uYb3E9pOvcpMjWcRAREaUTt0cyn46IyOEyJqirdHsTwa8cWGzvQIiIiNKIyy2ZT0dE5HAZE9S53N61ulxWvyQiIvJxezzI5rGRiMjRMiao01fqmDdARETk582pY1BHRORkGRPhVHKljoiIKAhz6oiInC9jgjqX1lw1Nztj/mQiIqKwWP2SiMj5MuZb/GC5GwBQLS/b5pEQERGlD67UERE5X8YEdUdc3qCuIJdBHRERkY45dUREzpcxQV2Fy7v9Mj8nY/5kIiKisNweD1fqiIgcLmMinHIGdURElCaEEG8LIbYLIRYr1z0ghNgkhJiv/XdyKsbicnP7JRGR02VMhKOv1OVlc/slERHZ7l0AI02uf0ZK2VP7b1wqBuL2SOSwMjQRkaNlXlDHlToiIrKZlPI3ALvsHgfgzanLZvVLIiJHy5hv8XKtUAq3XxIRURr7hxBiobY9s47ZHYQQ1woh5ggh5pSVlcX9hG4WSiEicryMiXC4UkdERGnuFQBtAfQEsAXAU2Z3klK+LqUskVKWNGjQIO4ndbFQChGR42VMhMNCKURElM6klNuklG4ppQfAGwD6puJ5uVJHROR8GRPhVLi8M5E52RnzJxMRkYMIIZooP54JYLHVfRPJxebjRESOl2P3AFJl7Y6DcHuk3cMgIiKCEOJjAMcDqC+E2AjgfgDHCyF6ApAASgFcl4qxcKWOiMj5MiaoG7toi91DICIiAgBIKS80ufqtlA8Eep867mIhInIyfosTERFlMK7UERE5H4M6IiKiDObyeJDN5uNERI6WMUFdr5a1Mbh9fbuHQURElFbWlB3E2IVMUSAicrKMyambt36P3UMgIiIiIiJKuIxZqSMiIiJzNfIzZo6XiKhKSnlQJ4RoIYSYLIRYKoRYIoS4JdVjICIiIq+a+Tk4r6SF3cMgIqI42DE15wLwbynln0KImgDmCiEmSCmXJusJpfT2p7uoX8tkPQUREZEjuaUEi18SETlbylfqpJRbpJR/apf3A1gGoFkyn7PS7Q3qmtWulsynISIichyPlMhmVEdE5Gi25tQJIYoB9AIwy+S2a4UQc4QQc8rKyuJ6niMuNwAgP4cphERERCqPBxCCQR0RkZPZFuUIIWoA+BLAP6WU+4y3Sylfl1KWSClLGjRoENdzlVd6ADCoIyIiMvJw+yURkePZEuUIIXLhDeg+klJ+leznK9dX6nKzk/1UREREjsLtl0REzmdH9UsB4C0Ay6SUT6fiOctdXKkjIiIyklLCI7n9kojI6eyIcgYCuBTAUCHEfO2/k5P5hAfLXQAY1BEREam04tDIZlBHRORoKW9pIKWcBiClR4/Hf1wBAFi8aR9GdmuSyqcmIiJKW24tquPuSyIiZ8uIpat+resCAE7r2dTmkRAREaUPjx7UMaojInK0jAjqqud7FyQb1sy3eSRERETpw+NNOUcWt18SETlaRgR1hyu91S8LWP2SiIjIR1+py86IswEioqorI77GyyvdEIKFUoiIiFS+7ZdcqSMicrSMiHIOV7pRkJPNks1EREQKffslj49ERM6WEUHdkUoPCnIz4k8lIiKKmG/7JWM6IiJHy4hI53Clm/l0REREBm5WvyQiqhJS3qfODu0a1vDNRhIREZFXblYWBrevj6a1qtk9FCIiikNGBHV/P66t3UMgIiJKO7UKc/HB1f3sHgYREcUpI7ZfEhERERERVVUM6oiIiIiIiByMQR0REREREZGDMagjIiIiIiJyMAZ1REREREREDsagjoiIiIiIyMEY1BERERERETkYgzoiIiIiIiIHY1BHRERERETkYAzqiIiIiIiIHIxBHRERERERkYMxqCMiIiIiInIwBnVEREREREQOJqSUdo8hLCFEGYB1cT5MfQA7EjCcqoSvSTC+Jub4ugTjaxIsEa9JKyllg0QMJhMk6PgI8P1shq9JML4mwfiaBONrYi6px0hHBHWJIISYI6UssXsc6YSvSTC+Jub4ugTjaxKMr4lz8d8uGF+TYHxNgvE1CcbXxFyyXxduvyQiIiIiInIwBnVEREREREQOlklB3et2DyAN8TUJxtfEHF+XYHxNgvE1cS7+2wXjaxKMr0kwvibB+JqYS+rrkjE5dURERERERFVRJq3UERERERERVTkZEdQJIUYKIVYIIVYLIUbbPZ5kEUK0EEJMFkIsFUIsEULcol1fVwgxQQixSvt/He16IYR4XntdFgoheiuPdbl2/1VCiMvt+psSRQiRLYSYJ4T4Qfu5tRBilva3fyqEyNOuz9d+Xq3dXqw8xp3a9SuEECNs+lMSRghRWwjxhRBiuRBimRBiQKa/V4QQt2qfncVCiI+FEAWZ9l4RQrwthNguhFisXJew94UQoo8QYpH2O88LIURq/0JSZcrxEeAxMhQeIwPx+GiOx8g0P0ZKKav0fwCyAawB0AZAHoAFALrYPa4k/a1NAPTWLtcEsBJAFwCPAxitXT8awGPa5ZMBjAcgAPQHMEu7vi6Av7T/19Eu17H774vztfkXgP8B+EH7+TMAF2iXXwVwvXb5BgCvapcvAPCpdrmL9t7JB9Bae09l2/13xfmavAfgGu1yHoDamfxeAdAMwFoA1ZT3yBWZ9l4BcCyA3gAWK9cl7H0BYLZ2X6H97kl2/82Z+h8y6Pio/b08Rlq/NjxGBr4ePD4GvyY8Rsr0PkZmwkpdXwCrpZR/SSkrAHwC4HSbx5QUUsotUso/tcv7ASyD90N4OrxfUND+f4Z2+XQA70uvmQBqCyGaABgBYIKUcpeUcjeACQBGpu4vSSwhRHMAowC8qf0sAAwF8IV2F+Nror9WXwAYpt3/dACfSCnLpZRrAayG973lSEKIWvB+Mb0FAFLKCinlHmT4ewVADoBqQogcAIUAtiDD3itSyt8A7DJcnZD3hXZbkZRypvQevd5XHotSL2OOjwCPkVZ4jAzE42NIPEam8TEyE4K6ZgA2KD9v1K6r0rRl7l4AZgFoJKXcot20FUAj7bLVa1PVXrNnAdwOwKP9XA/AHimlS/tZ/ft8f7t2+17t/lXtNWkNoAzAO9qWmzeFENWRwe8VKeUmAE8CWA/vgWovgLngewVI3PuimXbZeD3Zoyq+VyPCY2SAZ8FjpIrHRxM8RoaUFsfITAjqMo4QogaALwH8U0q5T71Ni/wzpuSpEOIUANullHPtHkuayYF3+8ArUspeAA7Cu2XAJwPfK3XgnVVrDaApgOpw/qxqwmXa+4KqHh4j/XiMNMXjowkeIyNj53sjE4K6TQBaKD83166rkoQQufAerD6SUn6lXb1NW9KF9v/t2vVWr01Ves0GAjhNCFEK79aioQCeg3cJPEe7j/r3+f527fZaAHaiar0mgHf2Z6OUcpb28xfwHsQy+b0yHMBaKWWZlLISwFfwvn8y/b0CJO59sUm7bLye7FEV36sh8RgZhMfIYDw+muMx0lpaHCMzIaj7A0B7rTpPHrzJmt/ZPKak0PYqvwVgmZTyaeWm7wDolXUuB/Ctcv1lWnWe/gD2asvHPwE4UQhRR5uZOVG7znGklHdKKZtLKYvh/befJKW8GMBkAOdodzO+JvprdY52f6ldf4FWzak1gPbwJrM6kpRyK4ANQoiO2lXDACxFBr9X4N1S0l8IUah9lvTXJKPfK5qEvC+02/YJIfprr/FlymNR6mXM8RHgMdIMj5HBeHy0xGOktfQ4Rso0qCST7P/grT6zEt4KO3fbPZ4k/p2D4F3yXQhgvvbfyfDuYf4FwCoAEwHU1e4vALykvS6LAJQoj3UVvMmrqwFcaffflqDX53j4K3u1gfdLZDWAzwHka9cXaD+v1m5vo/z+3dprtQJVoGIfgJ4A5mjvl2/grcCU0e8VAA8CWA5gMYAP4K3OlVHvFQAfw5svUQnvjPXViXxfACjRXt81AF4EIOz+mzP5v0w5Pmp/K4+RoV8fHiP9fwuPj+avC4+RaXyMFNoDEBERERERkQNlwvZLIiIiIiKiKotBHRERERERkYMxqCMiIiIiInIwBnVEREREREQOxqCOiIiIiIjIwRjUESWIEOJuIcQSIcRCIcR8IUS/JD7Xr0KIkmQ9PhERUSLxGEmUXDnh70JE4QghBgA4BUBvKWW5EKI+gDybh0VERGQ7HiOJko8rdUSJ0QTADillOQBIKXdIKTcLIe4TQvwhhFgshHhdCCEA3yziM0KIOUKIZUKIo4UQXwkhVgkhHtbuUyyEWC6E+Ei7zxdCiELjEwshThRC/C6E+FMI8bkQooZ2/RghxFJtVvTJFL4WREREKh4jiZKMQR1RYvwMoIUQYqUQ4mUhxHHa9S9KKY+WUnYDUA3emUpdhZSyBMCrAL4FcCOAbgCuEELU0+7TEcDLUsrOAPYBuEF9Um228x4Aw6WUvQHMAfAv7ffPBNBVStkdwMNJ+JuJiIgiwWMkUZIxqCNKACnlAQB9AFwLoAzAp0KIKwAMEULMEkIsAjAUQFfl177T/r8IwBIp5RZtFvMvAC202zZIKadrlz8EMMjw1P0BdAEwXQgxH8DlAFoB2AvgCIC3hBBnATiUqL+ViIgoGjxGEiUfc+qIEkRK6QbwK4BftQPUdQC6AyiRUm4QQjwAoED5lXLt/x7lsv6z/tmUxqcx/CwATJBSXmgcjxCiL4BhAM4B8A94D5hEREQpx2MkUXJxpY4oAYQQHYUQ7ZWregJYoV3eoe3hPyeGh26pJZgDwEUAphlunwlgoBCinTaO6kKIDtrz1ZJSjgNwK4AeMTw3ERFR3HiMJEo+rtQRJUYNAC8IIWoDcAFYDe82kz0AFgPYCuCPGB53BYAbhRBvA1gK4BX1RillmbaF5WMhRL529T0A9gP4VghRAO9M5b9ieG4iIqJE4DGSKMmElMaVaiJKB0KIYgA/aAnkREREpOExkigQt18SERERERE5GFfqiIiIiIiIHIwrdURERERERA7GoI6IiIiIiMjBGNQRERERERE5GIM6IiIiIiIiB2NQR0RERERE5GAM6oiIiIiIiBzs/wFJZIV3wGyrgwAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 1080x432 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Multiple chains (VMAP)"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "source": [
        "num_chains = 40\n",
        "initial_positions = {\"loc\": np.ones(num_chains), \"scale\": 2.0 * np.ones(num_chains)}\n",
        "initial_states = jax.vmap(hmc.new_state, in_axes=(0, None))(\n",
        "    initial_positions, potential\n",
        ")"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "source": [
        "\n",
        "\n",
        "def inference_loop_multiple_chains(rng_key, kernel, initial_state, num_samples, num_chains):\n",
        "    kernel_vmap = jax.vmap(kernel)\n",
        "    def one_step(states, rng_key):\n",
        "        keys = jax.random.split(rng_key, num_chains)\n",
        "        states, _ = kernel_vmap(keys, states)\n",
        "        return states, states\n",
        "\n",
        "    keys = jax.random.split(rng_key, num_samples)\n",
        "    _, states = jax.lax.scan(one_step, initial_state, keys)\n",
        "\n",
        "    return states"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "source": [
        "\n",
        "%%time\n",
        "\n",
        "states = inference_loop_multiple_chains(rng_key, hmc_kernel, initial_states, 2_000, num_chains)\n",
        "states.position[\"loc\"].block_until_ready()\n",
        "\n",
        "print(states.position[\"loc\"].shape)"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(2000, 40)\n",
            "CPU times: user 1.14 s, sys: 1.11 ms, total: 1.14 s\n",
            "Wall time: 1.82 s\n"
          ]
        }
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Multiple chains PMAP"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "source": [
        "num_devices = len(jax.local_devices())\n",
        "num_chains = 40\n",
        "initial_positions = {\"loc\": np.ones(num_devices),\n",
        "                    \"scale\": 2.0 * np.ones(num_devices)}                   \n",
        "initial_states = jax.pmap(hmc.new_state, in_axes=(0, None))(\n",
        "    initial_positions, potential\n",
        ")"
      ],
      "outputs": [
        {
          "output_type": "error",
          "ename": "TypeError",
          "evalue": "Argument '<function <lambda> at 0x7fd020503b80>' of type <class 'function'> is not a valid JAX type.",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
            "\u001b[0;32m/tmp/ipykernel_4242/256727136.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m initial_positions = {\"loc\": np.ones(num_devices),\n\u001b[1;32m      4\u001b[0m                     \"scale\": 2.0 * np.ones(num_devices)}                   \n\u001b[0;32m----> 5\u001b[0;31m initial_states = jax.pmap(hmc.new_state, in_axes=(0, None))(\n\u001b[0m\u001b[1;32m      6\u001b[0m     \u001b[0minitial_positions\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpotential\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m )\n",
            "    \u001b[0;31m[... skipping hidden 3 frame]\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.8/dist-packages/jax/_src/api.py\u001b[0m in \u001b[0;36m_check_arg\u001b[0;34m(arg)\u001b[0m\n\u001b[1;32m   2527\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_check_arg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2528\u001b[0m   \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTracer\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0m_valid_jaxtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2529\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"Argument '{arg}' of type {type(arg)} is not a valid JAX type.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2530\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2531\u001b[0m \u001b[0;31m# TODO(necula): this duplicates code in core.valid_jaxtype\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mTypeError\u001b[0m: Argument '<function <lambda> at 0x7fd020503b80>' of type <class 'function'> is not a valid JAX type."
          ]
        }
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "source": [
        "def inference_loop_pmap(rng_key, kernel, initial_state, num_samples, num_devices):\n",
        "    kernel_pmap = jax.pmap(kernel, axis_name='devices')\n",
        "    #kernel_vpmap = jax.vmap(kernel_pmap, axis_name='chains')\n",
        "    def one_step(states, rng_key):\n",
        "        keys = jax.random.split(rng_key, num_devices)\n",
        "        states, _ = kernel_pmap(keys, states)\n",
        "        return states, states\n",
        "\n",
        "    keys = jax.random.split(rng_key, num_samples)\n",
        "    _, states = jax.lax.scan(one_step, initial_state, keys)\n",
        "\n",
        "    return states"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "%%time\n",
        "\n",
        "states = inference_loop_pmap(rng_key, hmc_kernel, initial_states, 2_000, num_chains)\n",
        "states.position[\"loc\"].block_until_ready()\n",
        "\n",
        "print(states.position[\"loc\"].shape)"
      ],
      "outputs": [],
      "metadata": {}
    }
  ]
}